MATLAB等几何分析程序：NURBS基函数与刚度矩阵组装实战

简介这份资源是面向初学者的MATLAB等几何分析IGA入门程序包nliga_v1.3适合希望理解IGA原理、掌握NURBS数学概念并在MATLAB中实现数值模拟的学生与工程研究人员。压缩包共240个文件约12.34MB以194个.m脚本为主体配合.c/.cc/.h底层源码、.mexw64编译模块、.msh网格数据及makefile等构建文件覆盖NURBS基函数生成、模型构建、求解与后处理等环节。已有1415人学习下载。包内基础函数库可帮助读者理解控制点、权重与B样条曲面的构造方式求解器与示例案例则便于对照运行、观察位移与应力结果是衔接CAD几何与有限元分析的实用练习材料。1. 一份能跑通的 MATLAB 等几何分析程序到底解决了谁的痛点如果你做过结构仿真大概率经历过这种割裂CAD 里是 NURBS 曲面网格划分器里是线性四面体两边几何永远对不上光是把模型倒进求解器就得折腾半天倒进去还得担心网格畸变。等几何分析IGA的思路很直接——既然 CAD 用 NURBS 描述几何那就让 NURBS 基函数直接当形函数几何和离散统一省掉中间那层网格转换。这份 MATLAB 编制的 IGA 程序核心价值就在这儿它把 NURBS 基函数构造、控制点映射、刚度矩阵组装、边界条件施加这条链路用 MATLAB 脚本串起来了适合做板壳、二维弹性体这类中小规模问题的教学复现和算法验证。MATLAB 的好处是矩阵运算天然贴合有限元组装调试时能直接看基函数图像比 C 改一行编译半天友好得多。适合谁做结构仿真想入门 IGA 的研究生、需要快速验证新基函数或新积分方案的算法工程师以及手头有 NURBS 几何但不想碰商业软件二次开发的人。2. NURBS 基函数与等参映射程序的地基怎么搭2.1 从节点矢量到基函数值程序里最不能省的一步IGA 和传统有限元最大的差别就在基函数不是插值型的拉格朗日多项式而是 NURBS 基函数。程序里第一步永远是给定节点矢量knot vector比如一个 p 次、n 个控制点的一维情况节点矢量长度是 np1。MATLAB 里常见做法是先写一个openKnotVector函数生成 clamped 节点矢量再写basisFuns用 Cox-de Boor 递推算基函数值。这段代码是整个程序的根基函数算错后面刚度矩阵全是废的。function [N, dN] basisFuns(i, u, p, U) % i: 当前span索引, u: 参数点, p: 次数, U: 节点矢量 % 返回 p1 个非零基函数值 N 及其一阶导 dN N zeros(1, p1); dN zeros(1, p1); N(1) 1; for j 1:p left(j1) u - U(i1-j); right(j1) U(ij1) - u; saved 0; for r 0:j-1 denom right(r2) left(j-r1); temp N(r1) / denom; N(r1) saved right(r2) * temp; saved left(j-r1) * temp; end N(j1) saved; end % 导数部分略按标准递推补全 end这段是 Cox-de Boor 的标准实现left和right两个临时数组存的是当前参数点到左右节点的距离denom是分母saved用来滚动保存上一轮的值。参数说明i必须由findSpan函数先算出来不能随手给U必须是单调不减的重复度不能超过 p1否则分母为零直接 NaN。我一般会在调用前加一句assert(all(diff(U)0))做防御血泪经验是节点矢量手写错一位基函数图像看着还挺正常但积分结果偏 3% 以上查半天查不出来。2.2 等参映射与雅可比二维问题里最容易翻车的地方一维基函数算完二维就要做张量积。程序里通常用两个方向的节点矢量U、V和次数p、q控制点排成(n1)×(m1)的网格。每个高斯积分点上先分别算Nu、Nv再N kron(Nv, Nu)得到单元基函数向量。几何映射是x(ξ,η) Σ N_A(ξ,η) * P_A其中P_A是控制点坐标。雅可比矩阵J [∂x/∂ξ ∂x/∂η; ∂y/∂ξ ∂y/∂η]程序里必须用det(J)做积分权重同时inv(J)把物理坐标导数转成参数坐标导数。% 单元循环内高斯点 (xi, eta) 处 [Ru, dRu] basisFuns(spanU, xi, p, U); [Rv, dRv] basisFuns(spanV, eta, q, V); R kron(Rv, Ru); % 1 x (p1)(q1) dR [kron(Rv, dRu); kron(dRv, Ru)]; % 2 x nloc J dR * Pcoords; % 2 x 2, Pcoords 为单元控制点坐标 detJ det(J); if detJ 0 error(雅可比非正检查控制点顺序或节点矢量); end dRdx J \ dR; % 物理坐标导数这里Pcoords是当前单元涉及的控制点坐标矩阵行数等于局部基函数个数。detJ 0的判断是必须的控制点顺序反了或者节点矢量重复度给错雅可比会变负积分出来的刚度矩阵不对称求解器可能不报错但结果完全不对。常见做法是在单元循环里把detJ的最小值记下来跑完打印一遍如果接近零说明几何有退化得回头查控制点。另一个坑是kron的顺序kron(Rv, Ru)和kron(Ru, Rv)得到的基函数排列不一样必须和后面组装时的自由度编号严格对应否则刚度矩阵看着对称实际是错位的。2.3 刚度矩阵组装稀疏存储和自由度映射IGA 的组装逻辑和有限元几乎一样差别在单元连接关系不是网格给的而是由节点矢量的 span 推出来的。程序里一般先建一个elementConn矩阵每行存一个单元的控制点全局编号再建dofMap把控制点编号映射到自由度编号二维弹性问题每个控制点两个自由度。组装时用三元组(I, J, V)存稀疏矩阵最后sparse(I,J,V,nDof,nDof)一次性生成。nloc (p1)*(q1); Ke zeros(2*nloc, 2*nloc); for g 1:nGauss % ... 算 R, dRdx, detJ, 权重 w ... B zeros(3, 2*nloc); for a 1:nloc B(:, 2*a-1:2*a) [dRdx(1,a) 0; 0 dRdx(2,a); dRdx(2,a) dRdx(1,a)]; end Ke Ke B * D * B * detJ * w; end % 全局组装 rows dofMap(conn); cols rows; I [I; repmat(rows, 2*nloc, 1)]; ...D是材料矩阵平面应力问题就是那个含弹性模量和泊松比的 3×3 矩阵。B矩阵的构造顺序必须和dofMap一致常见错误是dRdx第一行对应 ξ 方向导数、第二行对应 η 方向但B里把 x/y 搞反结果位移场方向整体旋转 90 度。组装完建议先检查norm(K - K)IGA 刚度矩阵理论上对称不对称说明组装索引有重复或遗漏。稀疏矩阵的nDof要按控制点总数乘自由度算不是单元数这一点和有限元不同新手容易按单元节点数去开数组直接越界。3. 边界条件与求解控制点自由度怎么约束才不出错3.1 Dirichlet 边界条件的施加强施加和罚函数的选择IGA 里施加位移边界条件比有限元麻烦因为 NURBS 基函数不是插值型的控制点上的位移值不等于物理边界上的位移值。程序里常见两种做法一是强施加直接找到边界上对应的控制点把对应自由度从方程组里消掉二是罚函数法在刚度矩阵对角上加一个大数。这份程序大概率用的是强施加因为教学代码追求结果干净。强施加的步骤先根据参数域边界比如 ξ0 那条边反算出对应的控制点索引把这些控制点的自由度编号收集起来然后对刚度矩阵做行列置零、对角置一、右端项置位移值。MATLAB 里可以用K(bcDof,:) 0; K(:,bcDof) 0; K(bcDof,bcDof) speye(length(bcDof)); F(bcDof) u_bc;。注意K是稀疏矩阵时K(bcDof,:) 0这种操作效率不高控制点多了会明显变慢常见做法是先转成sparse再操作或者用逻辑索引一次性赋值。提示强施加时边界控制点的位移值不是随便给的如果边界是曲线控制点不一定落在物理边界上直接给常数值会引入误差。教学代码里通常假设边界是直线或者控制点恰好落在边界上实际工程问题要改用 Nitsche 方法或罚函数。3.2 高斯积分阶数与求解器选择IGA 单元内基函数是 p 次多项式刚度矩阵被积函数次数是 2p所以每个方向高斯点数至少取 p1 才能精确积分。程序里常见做法是nGauss p1但如果是变刚度或者几何非线性得往上加。我一般会跑两组积分阶数对比结果差超过 1e-6 就说明积分不够。求解器方面MATLAB 反斜杠对稀疏对称正定矩阵会自动选 Cholesky但 IGA 刚度矩阵条件数随次数升高恶化很快p 大于 4 时建议显式用ichol加预条件或者直接pcg带对角预条件。% 求解并检查残差 K (K K)/2; % 强制对称消除组装舍入误差 [L, U] lu(K); % 小规模直接分解 u U \ (L \ F); res norm(K*u - F) / norm(F); fprintf(相对残差: %.3e\n, res);(KK)/2这一步不是多余组装过程中浮点累加会让矩阵有 1e-15 量级的不对称对称求解器可能报错。res超过 1e-10 就要警惕要么是边界条件施加时行列置零破坏了对称性要么是某个单元雅可比接近零导致矩阵病态。常见误用是直接u K\F不检查残差跑出来位移场看着合理但应力云图在边界附近有异常振荡回头查半天才发现是某个控制点自由度没约束住。3.3 后处理从控制点位移到物理场解出来的是控制点位移要画云图得在参数域上取密集采样点用基函数插值回物理域。程序里一般写一个postProcess函数在 ξ、η 方向各取 50 个点算每个点的位移和应力再用scatter或trisurf画出来。应力计算要注意IGA 的应力在单元边界上连续但程序里如果按单元循环算边界点会被算两次画云图时会出现一条缝。常见做法是采样点避开单元边界或者用accumarray做平均。nSample 50; xiVec linspace(0, 1, nSample); [Xi, Eta] meshgrid(xiVec, xiVec); Ux zeros(size(Xi)); Uy zeros(size(Xi)); for a 1:numel(Xi) % 找 span, 算 R, 插值位移 [~, spanU] findSpan(...); R ...; Ux(a) R * uctrl(dofU); Uy(a) R * uctrl(dofV); end trisurf(delaunay(Xi(:), Eta(:)), Ux(:));采样密度nSample取 50 是经验值p 高的时候要加到 80 以上否则云图看着有棱角。delaunay在规则网格上效率不高可以直接用surf加参数网格但物理域不是矩形时得先映射。后处理这块最容易出的问题是位移插值用了错误的 span参数点落在单元边界上时findSpan返回的索引可能和相邻单元不一致导致云图在单元交界处跳变解决办法是采样点稍微偏离边界或者统一用左闭右开规则。4. 避坑与排查跑这份 IGA 程序时最常见的五类翻车4.1 基函数图像正常但结果偏几个百分点现象跑一个悬臂梁算例位移和理论解差 3% 到 5%基函数画出来光滑漂亮看不出问题。原因多半是节点矢量重复度给错clamped 节点矢量首尾重复度应该是 p1中间节点重复度不超过 p如果首尾只重复了 p 次边界基函数不满足 Kronecker 性质强施加边界条件时约束的不是真正的边界值。解决打印U的前 p1 个和后 p1 个元素确认首尾值相等且重复次数正确再检查findSpan在 u0 和 u1 处的返回值是否落在合法范围。4.2 刚度矩阵不对称或求解报奇异现象K组装完norm(K-K)是 1e-3 量级或者lu分解报矩阵奇异。原因通常是单元连接关系里控制点编号重复或越界kron顺序和dofMap不匹配或者某个单元detJ为负。解决先跑any(detJVec 0)定位坏单元再检查elementConn每行是否有重复编号最后核对dofMap的映射规则和B矩阵构造顺序是否一致。我一般会在组装前加一句assert(numel(unique(conn)) numel(conn))。4.3 边界条件施加后残差很大现象res在 1e-4 以上位移场在约束点附近有明显异常。原因可能是强施加时只置了行没置列或者置零操作把稀疏矩阵结构破坏了。解决确认K(bcDof,:) 0和K(:,bcDof) 0都执行了对角置一用speye而不是eye右端项在置零之后赋值。如果用的是罚函数罚因子取1e8 * max(abs(diag(K)))量级太小约束不住太大条件数爆炸。4.4 后处理云图在单元边界有缝现象位移云图在单元交界处出现一条细线或跳变。原因是采样点恰好落在单元边界findSpan返回的索引和相邻单元不一致插值结果不连续。解决采样点用linspace(0,1,nSample2)然后去掉首尾或者把参数点往单元内部偏移 1e-8。更稳妥的做法是按单元循环采样每个单元内部取点边界点用相邻单元平均值。4.5 高次基函数下求解慢到不可接受现象p 从 2 加到 4自由度没涨多少求解时间翻十倍。原因是刚度矩阵条件数随 p 指数恶化MATLAB 反斜杠的默认 pivoting 在病态矩阵上效率低。解决显式用ichol做不完全 Cholesky 预条件配合pcg容差设 1e-10或者对刚度矩阵做对角缩放K D\K/D把条件数压下来再求解。p 大于 5 时建议改用gmres带重启别硬扛。5. 进阶用法把这份程序改成自己的算法试验台跑通标准算例只是起点这份程序真正的价值在于它结构清晰改起来方便。我一般会做三件事第一把basisFuns抽成独立函数文件换基函数比如 T-spline 或 B 样条局部加密时只改这一个文件主程序不动第二在单元循环里加一个integrandHook函数句柄默认是标准B*D*B想试新积分方案或者加材料非线性时直接替换句柄不用改组装逻辑第三把postProcess里的采样和插值拆开采样点单独存成samplePts.mat换后处理方式时不用重算位移场。% 可替换积分被积函数的钩子写法 integrandHook (B, D, detJ, w) B * D * B * detJ * w; % 单元循环内 Ke Ke integrandHook(B, D, detJ, w); % 想试罚函数或者 Nitsche 时换一个句柄即可验证改动的正确性我习惯用 patch test取一个常应变场看数值解是否精确复现。IGA 的 patch test 比有限元严格因为基函数不是插值的常应变场下控制点位移得按基函数值反算不能直接给。具体做法是在边界上给线性位移场内部控制点位移用最小二乘拟合基函数值跑完看norm(u - u_exact)/norm(u_exact)是否在 1e-12 量级。这个测试能抓出基函数归一化错误、雅可比符号错误、组装索引错位这几类问题比对着云图肉眼看好使。另一个实用技巧是收敛性研究脚本化。写一个外层循环p 从 2 到 5每个 p 下把参数域均匀加密h-refinement记录自由度和能量范数误差画 log-log 图看收敛斜率。IGA 的 k-refinement先升次再加密和 h-refinement 收敛率不一样脚本化之后对比很直观。我一般会把结果存成convStudy.mat下次改算法直接重跑对比不用从头搭。注意k-refinement 时节点矢量插入要保证重复度不超过 p升次操作后原来的节点重复度含义变了直接复用旧矢量会出问题。常见做法是升次后重新生成 clamped 矢量再按新次数插节点。从那以后我每次拿到一份 IGA 代码都先跑 patch test 再跑收敛性两个都过了才敢往上面加新东西。这份 MATLAB 程序的好处是规模不大改坏了容易回滚适合当算法试验台用。希望帮到你。本文还有配套的精品资源点击获取