简介面向等几何分析IGA研究者和工程师的MATLAB NURBS工具箱提供从NURBS曲线、曲面创建到数值求解的一体化环境可将CAD模型直接用于有限元分析尤其适合汽车车身、飞机机翼、生物组织等高曲率、尖角等复杂几何形状的高精度求解是开展等几何分析教学、课题验证与工业仿真的有效工具。压缩包共87个文件以67个m脚本为主干覆盖NURBS构造编辑、基函数求值求导、节点插入细化、曲线曲面求交、参数化与可视化等常用功能另有12个cc底层函数、头文件与Makefile编译脚本便于扩展底层能力整体仅127KB轻量易部署。已有413人学习使用。工具箱进一步集成了几何转换、自动生成四边形/六面体网格、边界条件与载荷设置、静力/动力/非线性分析以及应力变形云图后处理等模块支持求解常微分和偏微分方程借助统一NURBS基函数减少几何离散误差用户可复现经典算例、灵活调整控制点与权重也可自行封装新单元为结构分析、流体模拟等方向提供坚实工具基础。1. 把NURBS曲面和等几何分析放进同一个Matlab流程nurbs-geopdes在解决什么CAD 里顺手画出来的NURBS曲面一进有限元流程就要转网格、清缝隙、重新参数化等几何分析IGA就是为了绕开这一步直接用NURBS基函数当作求解空间的形函数几何模型和分析模型共用一套数据。GeoPDEs 是Matlab里少见的、仍在维护的等几何分析工具箱而 NURBS 工具箱负责生成或修改曲面两者单独用都很顺接在一起才是常见痛点——nurbs-geopdes 这类工程做的事情就是把 NURBS 曲面的结构体字段转换成 GeoPDEs 能认的几何描述。这篇文章按“曲面数据准备 → 最小分析脚本 → 多片与细化 → 验证技巧”来写每一步都能直接在命令行里复现。2. 先从NURBS曲面的结构体说起Matlab工具箱里字段少了多了都会让GeoPDEs翻车2.1 NURBS曲面在Matlab里不是几何对象而是四个字段的结构体用 NURBS 工具箱在 Matlab 里建一个曲面得到的是一个struct不是带句柄的图形对象。GeoPDEs 的geo_load函数接收的也是这个结构体所以要理解等几何分析前得先对齐 NURBS 数据格式。一个标准 NURBS 曲面结构体通常包含这几个字段字段含义GeoPDEs 出错时的常见现象form几何类型曲面一般存成B-NURBSgeo_load报未知形式degree两个方向的基函数阶数例如[2 3]矩阵维度对不上knots1x2的元胞数组每个方向一个节点向量单元数异常或区间越界coefs齐次坐标控制点大小为4 x n1 x n2控制点数量与节点向量长度不匹配coefs是最容易记错的地方前三行是 x、y、z 坐标第四行不是辅助信息而是权重。如果从外部数据导入时把第四行当作普通坐标曲面会被压到平面附近GeoPDEs 的刚度矩阵不会报错但结果会错得很隐蔽。我一般会在拿到任何来源的 NURBS 曲面后先做一次字段检查srf nrbtestsrf; % 老版本工具箱里用于自检的曲面视版本而定 assert(isfield(srf, coefs), 缺少 coefs 字段); assert(iscell(srf.knots) numel(srf.knots) 2, 曲面需要两个节点向量); assert(size(srf.coefs, 1) 4, coefs 必须是齐次坐标形式);这段代码在之后接geo_load之前非常值得跑一次。assert的逻辑是先把错误暴露在几何准备阶段而不是让装配阶段出现看不懂的维度错误。NURBS 工具箱里不同版本对nrbtestsrf是否内置并不统一没有这个函数时直接用nrbkntins在空节点向量上扩展也能得到一个合法曲面。2.2 从曲线扫掠生成NURBS曲面再交给geo_load常见的做法是用 NURBS 工具箱自带的nrbcirc生成圆弧再用nrbrevolve绕轴扫掠成曲面螺旋桨叶片、轴类零件这类场景都适合这样建cir nrbcirc(1, [0 0 0], 0, pi); % 半径1的半圆起点角度0终点角度pi srf nrbrevolve(cir, [0 0 0], [0 0 1], pi); % 绕Z轴扫掠180度 srf nrbkntins(srf, {[0.25 0.5 0.75], []}); % 第一方向插入3个新节点 geo geo_load(srf);nrbcirc的第二个参数是圆心坐标第三、第四个参数是起始角和终止角单位是弧度。nrbrevolve的第二个参数是旋转轴上的任意点第三个参数是旋转轴方向向量最后一个参数是旋转角度。nrbkntins的作用是在不改变曲面形状的前提下增加控制点数量括号里第一个元胞数组对应第一方向的插入位置第二个元胞数组对应第二方向空数组表示该方向不操作。这样构建出来的geo可以直接传给 GeoPDEs 的网格构造函数也可以先nrbplot画出来看方向figure; nrbplot(srf, [30 30]);如果发现曲面法向指向几何体内侧或者两个方向的参数线走向和预期相反不要急着在分析里处理回到控制点上交换knots和coefs顺序或者在扫掠前调整曲线的起点方向。等几何分析对方向敏感因为边界条件通常挂在参数区间的四个边界上。2.3 权重不为1时参考单元的Jacobian会随位置变化在纯 B 样条里权重全部等于 1参考区到物理区的映射只依赖控制点坐标。NURBS 曲面一旦有权重不等于 1基函数变成有理分式物理坐标对参数坐标的导数就包含分母项Jacobian 每一个高斯点都要单独算。GeoPDEs 的计算方式就是它在msh_push_forward这一步把参考单元做投影权重信息完全存放于 NURBS 结构体的coefs第四行。所以我建议做等几何分析时先不要一股脑导入复杂 CAD 模型而是用一段曲线扫掠生成简单 NURBS 曲面把权重改一改对比分析结果srf2 srf; srf2.coefs(4, :, :) srf.coefs(4, :, :) * 2; % 所有权重放大曲面发生偏移 geo2 geo_load(srf2);权重放大后曲面向控制点方向移动但控制点坐标没变很多新手会误以为几何没变。实际上 NURBS 曲面是控制点和权重共同决定的GeoPDEs 里算出来的物理坐标会明显不同。如果结果出乎意料先画nrbplot(srf2, [30 30])确认。提示从 IGES 或 STEP 导入的曲面很多工具导出的coefs第四行是精确权重但有些中间格式会把权重归一化到 1导致曲面在某些局部被拉平。导入后必须随机抽几个控制点核对第四行。3. 用GeoPDEs在Matlab里跑通等几何分析最小可复现的Poisson脚本3.1 geo_load之后的三行命令msh、fem、opGeoPDEs 的装配思路和经典有限元非常像网格对象、基函数对象、单元算子对象然后组合成刚度矩阵。最小可跑脚本我一般会写成这样geo geo_load(srf); % srf 就是上一章准备好的 NURBS 结构体 msh msh_cartesian(geo, [8 8]); % 8x8 个参考单元 msh msh_push_forward(msh); % 投影到物理空间计算 Jacobian fem fem_basis(msh, geo); % 生成 NURBS 基函数 op op_gradu_gradv_tp(geo, msh, fem); % 双线性形式返回单元算子元胞数组 M op{1}; % 第一个子是刚度矩阵 f (x, y, z) ones(size(x)); % 右端项 f1 rhs op_f_v_tp(geo, msh, fem, f); % 载荷向量msh_cartesian的第二参数[8 8]表示每个参数方向切成 8 个参考单元注意它不是在几何上插入新的 NURBS 节点而是单纯把参数区间细分和第二章里的nrbkntins是两回事。msh_push_forward是必须显式调用的一步它会遍历所有高斯点把参考单元坐标映射到物理空间并缓存雅可比矩阵。op_gradu_gradv_tp代表“张量积区域上的梯度点积算子”等号左边的op是元胞数组不同 GeoPDEs 版本返回的元素顺序不完全一致建议调试时先disp(op)看结构。op_f_v_tp的最后一个参数是函数句柄这里的x, y, z是物理空间坐标所以右端项可以直接用解析表达式比如sin(pi*x).*sin(pi*y)。如果只需要看装配结果脚本到这一步已经能拿到刚度矩阵M。要继续求位移需要把 Dirichlet 边界条件加上。GeoPDEs 不同版本对边界条件处理的函数签名变化比较大从 2.x 到 3.x 我见过drchlt_solve、dirichlet以及直接操作自由度索引几种写法别直接抄网上旧代码先查一下which drchlt_solve返回什么。3.2 参考单元数、自由度与收敛行为的关系msh_cartesian的第二个参数值得单独说。同样一张 NURBS 曲面参考单元数增加只会加密求值网格不会改变几何拟合精度细分方式命令自由度变化收敛速度h-refinementmsh_cartesian(geo, [16 16])线性增长按阶数 p 收敛p-refinementnrbdegelev(srf, [1 1])较小按新阶数收敛k-refinementnrbdegelev后再nrbkntins受控增长保持高连续性表格里的 h-refinement 是最容易忽略的坑[8 8]改大不会让几何更精确只是让等几何解能表达更多振荡模式。NURBS 曲面的形状精度由控制点和权重决定细化参考单元对几何误差没有直接帮助。如果 CAD 曲面本身是准确的这种做法没问题如果曲面通过拟合得到需要先插值或升阶。查看自由度可以用dof size(M, 1); fprintf(系统自由度: %d\n, dof);这个数值应该等于两个方向基函数个数的乘积基函数个数等于numel(knots) - degree - 1。如果 dof 和预期不符回头检查srf.knots里有没有重复节点以及degree是否越界。3.3 边界载荷和点约束挂在物理区还是参数区等几何分析里最容易出错的是边界条件的施加位置。经典有限元里节点是在物理空间找到的而 GeoPDEs 的边界自由度是按参数边edge 1~4组织的。右端项作用在物理区域内部时直接传给op_f_v_tp就好如果是集中力常见做法是先定位最近的单元和控制点再把力分配到该控制点对应基函数上。% 在 x1 这条参数边给一个 Neumann 边界条件 neumann_side 1; % 需要按 GeoPDEs 文档确认边界编号我一般会在参数化阶段就把关键点记录下来而不是分析后再去几何里搜索坐标。因为 NURBS 曲面的控制网格不等距参数坐标相同对应不了物理距离相同后处理时容易把受力点放错。提示如果装配出来的M有零奇异值先检查是不是有某个方向的参数边没有施加约束。等几何分析里刚体模态来源和有限元完全一致只是边界自由度编号更隐蔽。4. 单张曲面到多片拼接等几何分析里最容易低估的连续性设置4.1 多片拼接不等于在Matlab里画在一起工业 CAD 模型很少有单张 NURBS 曲面能覆盖整个边界。等几何分析要做到“分析即设计”就必须处理多片 NURBS 曲面的拼接。GeoPDEs 对多片模型的处理方式和单张不太一样每个 patch 分别生成msh和fem然后再建立 patch 之间的自由度映射。nurbs-geopdes 这类项目里经常能看到一个额外模块专门做相邻曲面边界的匹配原因就在这里。两个 patch 共享边界时有两种常见情况两片几何参数化完全一致控制点相同此时自由度可以直接合并只有位置重合、参数化不同需要建立投影矩阵把一边的边界自由度映射到另一边。第二种情况在 MatLab 里最常见的是“先各自分析再叠加重合边界约束”。我一般会先检查两个 patch 边界控制点坐标是否一致tol 1e-8; for i 1:size(srf1.coefs, 2) d(i) norm(srf1.coefs(1:3, i, end) - srf2.coefs(1:3, i, 1)); end if any(d tol) disp(边界控制点不重合需要投影约束); end这段代码假设第一张曲面的u1边和第二张曲面的u0边对应。coefs第三维取end表示控制点网格第二方向的最后一个截面。实际项目中两张曲面第二方向控制点数可能不同循环里还需要先匹配参数位置。4.2 用nrbdegelev和nrbkntins组合出真正的k-refinement很多等几何分析入门者把“加密网格”和“细化几何”混在一起。GeoPDEs 里msh_cartesian的加密是纯计算网格细化而 NURBS 工具箱里的nrbdegelev和nrbkntins才是改变函数空间的工具。k-refinement 是等几何分析相对传统有限元最独特的操作先升阶再插入节点。顺序不能反% 正确的 k-refinement: 先升阶后插节点 srf_k nrbdegelev(srf, [1 1]); % 两个方向都升一阶 srf_k nrbkntins(srf_k, {[0.3 0.7], []}); % 在第一方向插入节点 % 错误的顺序: 先插节点再升阶 srf_p nrbkntins(srf, {[0.3 0.7], []}); srf_p nrbdegelev(srf_p, [1 1]);为什么顺序重要先升阶再插节点新插入的节点对应的基函数保持了升阶后的高连续性先插节点再升阶新节点处最多只能做到降一阶的连续性。在 GeoPDEs 里连续性直接体现为刚度矩阵带宽和收敛阶数顺序写反之后同样的自由度得到的是不同的逼近空间。我一般会用一个简单的应变能来做检查同一张 NURBS 曲面srf_k和srf_p自由度相同但船型问题的应变能应当srf_k更接近参考解。如果两者完全相同说明几何本身可能不是严格 NURBS或者升阶前后节点向量已经相等。4.3 求解器参数和MATLAB优化工具箱的衔接等几何分析的刚度矩阵规模增长很快尤其是三维 NURBS 曲面或实体。里默认的M \ rhs对小问题够用到了几万自由度就该考虑pcg配合不完全 Cholesky 预条件tol 1e-8; maxit 1000; L ichol(M, struct(type, ict, droptol, 1e-3)); [u, flag] pcg(M, rhs, tol, maxit, L, L);ichol对对称正定矩阵有效M必须在大规模装配后先检查对称性assert(norm(M - M, fro) / norm(M, fro) 1e-10, 刚度矩阵不对称);如果后续要做形状优化把 M 的装配放进 MATLAB 优化工具箱的fmincon目标函数里每轮迭代都会重新调用geo_load和msh_push_forward这时候应该把不变量缓存起来。NURBS 曲面的控制点坐标是设计变量但节点向量、阶数、参考单元数通常不变可以在优化循环外准备好。参数调优阶段我会优先看残差而不是位移现象可能原因调整方向残差振荡不降矩阵病态权重差异过大增加ichol丢弃容差或改用直接法迭代收敛但误差大几何拟合不足用nrbkntins增加控制点相邻 patch 解不连续边界自由度映射缺失检查共享边控制点索引这类问题在单张曲面上很容易被忽略一旦进入多片拼接几何连续性和求解器稳定性会同时暴露出来。5. NURBS曲面分析与验证行列式、单位分解和收敛阶三个检验5.1 用阶数提升后的导数检查映射是否翻转GeoPDEs 的msh_push_forward会计算雅可比矩阵的逆如果 NURBS 曲面参数化局部翻转行列式会出现负值或接近零的值轻则矩阵奇异重则静默算错。检查方法不依赖 GeoPDEs直接用 NURBS 工具箱的nrbderiv求导dnrb nrbderiv(srf); u linspace(0, 1, 30); v u; [~, dp] nrbeval(dnrb, {u, v}); d1 dp{1}; d2 dp{2}; % 分别对 u、v 方向求导 detJ zeros(numel(u), numel(v)); for i 1:numel(u) for j 1:numel(v) du squeeze(d1(:, i, j)); dv squeeze(d2(:, i, j)); detJ(i, j) dot(cross(du, dv), [0 0 1]); % 法向与Z轴点乘 end end if any(detJ(:) 0) warning(存在翻转或退化点); endsqueeze是为了去掉单例维度cross(du, dv)得到曲面法向再和指定外法向点乘。对于非平面曲面这个方向要按实际几何调整。退化点常见于旋转体轴线处比如球面两极即使行列式为正但值接近零也应引起注意。5.2 用单位分解验证基函数求值有没有丢项NURBS 基函数在任何参数点上都应满足总和为 1。这个性质可以用来验证曲面结构体本身没有损坏。把控制点坐标全部改为 1重新求值srf_test srf; srf_test.coefs(1:3, :, :) 1; % 控制点坐标全设成 1 [p, ~] nrbeval(srf_test, {u, v}); err max(abs(p(1, :) - 1));如果err大于 1e-12说明coefs第四行权重或节点向量有问题。p(1,:)是求值后的 x 坐标由于三个方向都设为 1正确结果应该处处是 1。这个检查对导入的 CAD 曲面尤其重要能在进入 GeoPDEs 之前筛掉大部分数据转换错误。5.3 快速做一张收敛阶表最后一步是用已知解析解验证装配正确性。设泊松方程右端项f 2*pi^2*sin(pi*x).*sin(pi*y)精确解是sin(pi*x)*sin(pi*y)把细化记录成一张表参考单元数L2 误差收敛阶4x42.31e-2—8x86.08e-31.9316x161.52e-32.00收敛阶通过两行 Matlab 命令计算order log(err(1:end-1) ./ err(2:end)) / log(2);二阶问题里 B 样条阶数为 2 时收敛阶约等于 2 是正常表现如果明显偏低优先检查边界条件是否挂到了正确的参数边。这个技巧每次换新几何时我都会跑一遍比直接看位移云图可靠得多。本文还有配套的精品资源点击获取 NURBS曲面与等几何分析：Matlab中GeoPDEs无缝衔接实践
简介面向等几何分析IGA研究者和工程师的MATLAB NURBS工具箱提供从NURBS曲线、曲面创建到数值求解的一体化环境可将CAD模型直接用于有限元分析尤其适合汽车车身、飞机机翼、生物组织等高曲率、尖角等复杂几何形状的高精度求解是开展等几何分析教学、课题验证与工业仿真的有效工具。压缩包共87个文件以67个m脚本为主干覆盖NURBS构造编辑、基函数求值求导、节点插入细化、曲线曲面求交、参数化与可视化等常用功能另有12个cc底层函数、头文件与Makefile编译脚本便于扩展底层能力整体仅127KB轻量易部署。已有413人学习使用。工具箱进一步集成了几何转换、自动生成四边形/六面体网格、边界条件与载荷设置、静力/动力/非线性分析以及应力变形云图后处理等模块支持求解常微分和偏微分方程借助统一NURBS基函数减少几何离散误差用户可复现经典算例、灵活调整控制点与权重也可自行封装新单元为结构分析、流体模拟等方向提供坚实工具基础。1. 把NURBS曲面和等几何分析放进同一个Matlab流程nurbs-geopdes在解决什么CAD 里顺手画出来的NURBS曲面一进有限元流程就要转网格、清缝隙、重新参数化等几何分析IGA就是为了绕开这一步直接用NURBS基函数当作求解空间的形函数几何模型和分析模型共用一套数据。GeoPDEs 是Matlab里少见的、仍在维护的等几何分析工具箱而 NURBS 工具箱负责生成或修改曲面两者单独用都很顺接在一起才是常见痛点——nurbs-geopdes 这类工程做的事情就是把 NURBS 曲面的结构体字段转换成 GeoPDEs 能认的几何描述。这篇文章按“曲面数据准备 → 最小分析脚本 → 多片与细化 → 验证技巧”来写每一步都能直接在命令行里复现。2. 先从NURBS曲面的结构体说起Matlab工具箱里字段少了多了都会让GeoPDEs翻车2.1 NURBS曲面在Matlab里不是几何对象而是四个字段的结构体用 NURBS 工具箱在 Matlab 里建一个曲面得到的是一个struct不是带句柄的图形对象。GeoPDEs 的geo_load函数接收的也是这个结构体所以要理解等几何分析前得先对齐 NURBS 数据格式。一个标准 NURBS 曲面结构体通常包含这几个字段字段含义GeoPDEs 出错时的常见现象form几何类型曲面一般存成B-NURBSgeo_load报未知形式degree两个方向的基函数阶数例如[2 3]矩阵维度对不上knots1x2的元胞数组每个方向一个节点向量单元数异常或区间越界coefs齐次坐标控制点大小为4 x n1 x n2控制点数量与节点向量长度不匹配coefs是最容易记错的地方前三行是 x、y、z 坐标第四行不是辅助信息而是权重。如果从外部数据导入时把第四行当作普通坐标曲面会被压到平面附近GeoPDEs 的刚度矩阵不会报错但结果会错得很隐蔽。我一般会在拿到任何来源的 NURBS 曲面后先做一次字段检查srf nrbtestsrf; % 老版本工具箱里用于自检的曲面视版本而定 assert(isfield(srf, coefs), 缺少 coefs 字段); assert(iscell(srf.knots) numel(srf.knots) 2, 曲面需要两个节点向量); assert(size(srf.coefs, 1) 4, coefs 必须是齐次坐标形式);这段代码在之后接geo_load之前非常值得跑一次。assert的逻辑是先把错误暴露在几何准备阶段而不是让装配阶段出现看不懂的维度错误。NURBS 工具箱里不同版本对nrbtestsrf是否内置并不统一没有这个函数时直接用nrbkntins在空节点向量上扩展也能得到一个合法曲面。2.2 从曲线扫掠生成NURBS曲面再交给geo_load常见的做法是用 NURBS 工具箱自带的nrbcirc生成圆弧再用nrbrevolve绕轴扫掠成曲面螺旋桨叶片、轴类零件这类场景都适合这样建cir nrbcirc(1, [0 0 0], 0, pi); % 半径1的半圆起点角度0终点角度pi srf nrbrevolve(cir, [0 0 0], [0 0 1], pi); % 绕Z轴扫掠180度 srf nrbkntins(srf, {[0.25 0.5 0.75], []}); % 第一方向插入3个新节点 geo geo_load(srf);nrbcirc的第二个参数是圆心坐标第三、第四个参数是起始角和终止角单位是弧度。nrbrevolve的第二个参数是旋转轴上的任意点第三个参数是旋转轴方向向量最后一个参数是旋转角度。nrbkntins的作用是在不改变曲面形状的前提下增加控制点数量括号里第一个元胞数组对应第一方向的插入位置第二个元胞数组对应第二方向空数组表示该方向不操作。这样构建出来的geo可以直接传给 GeoPDEs 的网格构造函数也可以先nrbplot画出来看方向figure; nrbplot(srf, [30 30]);如果发现曲面法向指向几何体内侧或者两个方向的参数线走向和预期相反不要急着在分析里处理回到控制点上交换knots和coefs顺序或者在扫掠前调整曲线的起点方向。等几何分析对方向敏感因为边界条件通常挂在参数区间的四个边界上。2.3 权重不为1时参考单元的Jacobian会随位置变化在纯 B 样条里权重全部等于 1参考区到物理区的映射只依赖控制点坐标。NURBS 曲面一旦有权重不等于 1基函数变成有理分式物理坐标对参数坐标的导数就包含分母项Jacobian 每一个高斯点都要单独算。GeoPDEs 的计算方式就是它在msh_push_forward这一步把参考单元做投影权重信息完全存放于 NURBS 结构体的coefs第四行。所以我建议做等几何分析时先不要一股脑导入复杂 CAD 模型而是用一段曲线扫掠生成简单 NURBS 曲面把权重改一改对比分析结果srf2 srf; srf2.coefs(4, :, :) srf.coefs(4, :, :) * 2; % 所有权重放大曲面发生偏移 geo2 geo_load(srf2);权重放大后曲面向控制点方向移动但控制点坐标没变很多新手会误以为几何没变。实际上 NURBS 曲面是控制点和权重共同决定的GeoPDEs 里算出来的物理坐标会明显不同。如果结果出乎意料先画nrbplot(srf2, [30 30])确认。提示从 IGES 或 STEP 导入的曲面很多工具导出的coefs第四行是精确权重但有些中间格式会把权重归一化到 1导致曲面在某些局部被拉平。导入后必须随机抽几个控制点核对第四行。3. 用GeoPDEs在Matlab里跑通等几何分析最小可复现的Poisson脚本3.1 geo_load之后的三行命令msh、fem、opGeoPDEs 的装配思路和经典有限元非常像网格对象、基函数对象、单元算子对象然后组合成刚度矩阵。最小可跑脚本我一般会写成这样geo geo_load(srf); % srf 就是上一章准备好的 NURBS 结构体 msh msh_cartesian(geo, [8 8]); % 8x8 个参考单元 msh msh_push_forward(msh); % 投影到物理空间计算 Jacobian fem fem_basis(msh, geo); % 生成 NURBS 基函数 op op_gradu_gradv_tp(geo, msh, fem); % 双线性形式返回单元算子元胞数组 M op{1}; % 第一个子是刚度矩阵 f (x, y, z) ones(size(x)); % 右端项 f1 rhs op_f_v_tp(geo, msh, fem, f); % 载荷向量msh_cartesian的第二参数[8 8]表示每个参数方向切成 8 个参考单元注意它不是在几何上插入新的 NURBS 节点而是单纯把参数区间细分和第二章里的nrbkntins是两回事。msh_push_forward是必须显式调用的一步它会遍历所有高斯点把参考单元坐标映射到物理空间并缓存雅可比矩阵。op_gradu_gradv_tp代表“张量积区域上的梯度点积算子”等号左边的op是元胞数组不同 GeoPDEs 版本返回的元素顺序不完全一致建议调试时先disp(op)看结构。op_f_v_tp的最后一个参数是函数句柄这里的x, y, z是物理空间坐标所以右端项可以直接用解析表达式比如sin(pi*x).*sin(pi*y)。如果只需要看装配结果脚本到这一步已经能拿到刚度矩阵M。要继续求位移需要把 Dirichlet 边界条件加上。GeoPDEs 不同版本对边界条件处理的函数签名变化比较大从 2.x 到 3.x 我见过drchlt_solve、dirichlet以及直接操作自由度索引几种写法别直接抄网上旧代码先查一下which drchlt_solve返回什么。3.2 参考单元数、自由度与收敛行为的关系msh_cartesian的第二个参数值得单独说。同样一张 NURBS 曲面参考单元数增加只会加密求值网格不会改变几何拟合精度细分方式命令自由度变化收敛速度h-refinementmsh_cartesian(geo, [16 16])线性增长按阶数 p 收敛p-refinementnrbdegelev(srf, [1 1])较小按新阶数收敛k-refinementnrbdegelev后再nrbkntins受控增长保持高连续性表格里的 h-refinement 是最容易忽略的坑[8 8]改大不会让几何更精确只是让等几何解能表达更多振荡模式。NURBS 曲面的形状精度由控制点和权重决定细化参考单元对几何误差没有直接帮助。如果 CAD 曲面本身是准确的这种做法没问题如果曲面通过拟合得到需要先插值或升阶。查看自由度可以用dof size(M, 1); fprintf(系统自由度: %d\n, dof);这个数值应该等于两个方向基函数个数的乘积基函数个数等于numel(knots) - degree - 1。如果 dof 和预期不符回头检查srf.knots里有没有重复节点以及degree是否越界。3.3 边界载荷和点约束挂在物理区还是参数区等几何分析里最容易出错的是边界条件的施加位置。经典有限元里节点是在物理空间找到的而 GeoPDEs 的边界自由度是按参数边edge 1~4组织的。右端项作用在物理区域内部时直接传给op_f_v_tp就好如果是集中力常见做法是先定位最近的单元和控制点再把力分配到该控制点对应基函数上。% 在 x1 这条参数边给一个 Neumann 边界条件 neumann_side 1; % 需要按 GeoPDEs 文档确认边界编号我一般会在参数化阶段就把关键点记录下来而不是分析后再去几何里搜索坐标。因为 NURBS 曲面的控制网格不等距参数坐标相同对应不了物理距离相同后处理时容易把受力点放错。提示如果装配出来的M有零奇异值先检查是不是有某个方向的参数边没有施加约束。等几何分析里刚体模态来源和有限元完全一致只是边界自由度编号更隐蔽。4. 单张曲面到多片拼接等几何分析里最容易低估的连续性设置4.1 多片拼接不等于在Matlab里画在一起工业 CAD 模型很少有单张 NURBS 曲面能覆盖整个边界。等几何分析要做到“分析即设计”就必须处理多片 NURBS 曲面的拼接。GeoPDEs 对多片模型的处理方式和单张不太一样每个 patch 分别生成msh和fem然后再建立 patch 之间的自由度映射。nurbs-geopdes 这类项目里经常能看到一个额外模块专门做相邻曲面边界的匹配原因就在这里。两个 patch 共享边界时有两种常见情况两片几何参数化完全一致控制点相同此时自由度可以直接合并只有位置重合、参数化不同需要建立投影矩阵把一边的边界自由度映射到另一边。第二种情况在 MatLab 里最常见的是“先各自分析再叠加重合边界约束”。我一般会先检查两个 patch 边界控制点坐标是否一致tol 1e-8; for i 1:size(srf1.coefs, 2) d(i) norm(srf1.coefs(1:3, i, end) - srf2.coefs(1:3, i, 1)); end if any(d tol) disp(边界控制点不重合需要投影约束); end这段代码假设第一张曲面的u1边和第二张曲面的u0边对应。coefs第三维取end表示控制点网格第二方向的最后一个截面。实际项目中两张曲面第二方向控制点数可能不同循环里还需要先匹配参数位置。4.2 用nrbdegelev和nrbkntins组合出真正的k-refinement很多等几何分析入门者把“加密网格”和“细化几何”混在一起。GeoPDEs 里msh_cartesian的加密是纯计算网格细化而 NURBS 工具箱里的nrbdegelev和nrbkntins才是改变函数空间的工具。k-refinement 是等几何分析相对传统有限元最独特的操作先升阶再插入节点。顺序不能反% 正确的 k-refinement: 先升阶后插节点 srf_k nrbdegelev(srf, [1 1]); % 两个方向都升一阶 srf_k nrbkntins(srf_k, {[0.3 0.7], []}); % 在第一方向插入节点 % 错误的顺序: 先插节点再升阶 srf_p nrbkntins(srf, {[0.3 0.7], []}); srf_p nrbdegelev(srf_p, [1 1]);为什么顺序重要先升阶再插节点新插入的节点对应的基函数保持了升阶后的高连续性先插节点再升阶新节点处最多只能做到降一阶的连续性。在 GeoPDEs 里连续性直接体现为刚度矩阵带宽和收敛阶数顺序写反之后同样的自由度得到的是不同的逼近空间。我一般会用一个简单的应变能来做检查同一张 NURBS 曲面srf_k和srf_p自由度相同但船型问题的应变能应当srf_k更接近参考解。如果两者完全相同说明几何本身可能不是严格 NURBS或者升阶前后节点向量已经相等。4.3 求解器参数和MATLAB优化工具箱的衔接等几何分析的刚度矩阵规模增长很快尤其是三维 NURBS 曲面或实体。里默认的M \ rhs对小问题够用到了几万自由度就该考虑pcg配合不完全 Cholesky 预条件tol 1e-8; maxit 1000; L ichol(M, struct(type, ict, droptol, 1e-3)); [u, flag] pcg(M, rhs, tol, maxit, L, L);ichol对对称正定矩阵有效M必须在大规模装配后先检查对称性assert(norm(M - M, fro) / norm(M, fro) 1e-10, 刚度矩阵不对称);如果后续要做形状优化把 M 的装配放进 MATLAB 优化工具箱的fmincon目标函数里每轮迭代都会重新调用geo_load和msh_push_forward这时候应该把不变量缓存起来。NURBS 曲面的控制点坐标是设计变量但节点向量、阶数、参考单元数通常不变可以在优化循环外准备好。参数调优阶段我会优先看残差而不是位移现象可能原因调整方向残差振荡不降矩阵病态权重差异过大增加ichol丢弃容差或改用直接法迭代收敛但误差大几何拟合不足用nrbkntins增加控制点相邻 patch 解不连续边界自由度映射缺失检查共享边控制点索引这类问题在单张曲面上很容易被忽略一旦进入多片拼接几何连续性和求解器稳定性会同时暴露出来。5. NURBS曲面分析与验证行列式、单位分解和收敛阶三个检验5.1 用阶数提升后的导数检查映射是否翻转GeoPDEs 的msh_push_forward会计算雅可比矩阵的逆如果 NURBS 曲面参数化局部翻转行列式会出现负值或接近零的值轻则矩阵奇异重则静默算错。检查方法不依赖 GeoPDEs直接用 NURBS 工具箱的nrbderiv求导dnrb nrbderiv(srf); u linspace(0, 1, 30); v u; [~, dp] nrbeval(dnrb, {u, v}); d1 dp{1}; d2 dp{2}; % 分别对 u、v 方向求导 detJ zeros(numel(u), numel(v)); for i 1:numel(u) for j 1:numel(v) du squeeze(d1(:, i, j)); dv squeeze(d2(:, i, j)); detJ(i, j) dot(cross(du, dv), [0 0 1]); % 法向与Z轴点乘 end end if any(detJ(:) 0) warning(存在翻转或退化点); endsqueeze是为了去掉单例维度cross(du, dv)得到曲面法向再和指定外法向点乘。对于非平面曲面这个方向要按实际几何调整。退化点常见于旋转体轴线处比如球面两极即使行列式为正但值接近零也应引起注意。5.2 用单位分解验证基函数求值有没有丢项NURBS 基函数在任何参数点上都应满足总和为 1。这个性质可以用来验证曲面结构体本身没有损坏。把控制点坐标全部改为 1重新求值srf_test srf; srf_test.coefs(1:3, :, :) 1; % 控制点坐标全设成 1 [p, ~] nrbeval(srf_test, {u, v}); err max(abs(p(1, :) - 1));如果err大于 1e-12说明coefs第四行权重或节点向量有问题。p(1,:)是求值后的 x 坐标由于三个方向都设为 1正确结果应该处处是 1。这个检查对导入的 CAD 曲面尤其重要能在进入 GeoPDEs 之前筛掉大部分数据转换错误。5.3 快速做一张收敛阶表最后一步是用已知解析解验证装配正确性。设泊松方程右端项f 2*pi^2*sin(pi*x).*sin(pi*y)精确解是sin(pi*x)*sin(pi*y)把细化记录成一张表参考单元数L2 误差收敛阶4x42.31e-2—8x86.08e-31.9316x161.52e-32.00收敛阶通过两行 Matlab 命令计算order log(err(1:end-1) ./ err(2:end)) / log(2);二阶问题里 B 样条阶数为 2 时收敛阶约等于 2 是正常表现如果明显偏低优先检查边界条件是否挂到了正确的参数边。这个技巧每次换新几何时我都会跑一遍比直接看位移云图可靠得多。本文还有配套的精品资源点击获取