简介这份名为 nliga_v1.3 的 MATLAB 程序包是一套完整的等几何分析IGA入门实现框架面向计算力学初学者及希望将 CAD 建模与数值仿真衔接的研究人员。程序包以 NURBS 为基函数将几何建模与有限元分析统一在同一数学描述中可有效减少传统网格转换带来的几何误差。压缩包共 240 个文件总大小约 12.34MB主体为 194 个 m 脚本覆盖基函数生成、模型装配与后处理另有 12 个 cc、7 个 c 及 2 个 h 文件提供底层数值内核配合 11 个 msh 网格示例与 6 个 mexw64 已编译接口便于直接调用与二次开发。目前已有 1400 余人学习下载适合从零理解 IGA 原理并动手实践。借助内置示例与教程脚本使用者可掌握 NURBS 几何构建、边界条件设置、求解器调用及结果可视化形成一套可迁移的等几何分析工作流为后续深入研究结构力学或流体问题打下扎实基础。1. 等几何分析上手让结构仿真省掉整条网格生成流水线做结构仿真的人都有共识从CAD模型到有限元网格这一步常常比求解本身还费时间。模型复杂一点几何清理、网格加密、质量检查一套流程下来半天就没了。IGA等几何分析Isogeometric Analysis的思路是把这个问题从根本上绕过去不生成网格直接把CAD里的NURBS基函数当作分析形函数几何模型就是计算模型。这份资源是一套用MATLAB编制的IGA程序从NURBS基函数计算、几何映射、刚度矩阵组装到边界条件求解全流程覆盖跑一遍就能算出二维线弹性问题的位移场和应力场。适合有限元基础扎实、想搞懂等几何分析到底怎么写成代码的工程师也适合论文里需要参数化仿真复现的研究生。下面我按这套程序的核心模块往下拆重点说编码思路和可以复用的实现片段。2. 核心模块拆解NURBS基函数与几何映射的完整实现2.1 从递推公式到MATLAB代码B样条基函数等几何分析和传统有限元第一个分岔路口是形函数。传统有限元用拉格朗日多项式IGA改用B样条递推公式。B样条基函数的计算结构是一棵递归树p阶基函数由两个p-1阶基函数线性组合而来写进MATLAB非常紧凑。下面这段代码是这套程序里最基础的工具函数计算单条基函数在某个参数坐标处的值。function N bspline_basis(i, p, xi, Xi) % 计算第i个p阶B样条基函数在参数坐标xi处的值 % i: 基函数序号从1开始 % p: 基函数阶数p2表示二次样条 % xi: 参数坐标取值落在节点向量定义区间内 % Xi: 节点向量形如 [0 0 0 1 1 1] % 返回: 基函数值标量 if p 0 % 零阶基函数就是区间指示函数 if xi Xi(i) xi Xi(i1) N 1; else N 0; end return; end % 第一个非零项左侧基函数贡献 w1 0; den1 Xi(ip) - Xi(i); if den1 0 w1 (xi - Xi(i)) / den1 * bspline_basis(i, p-1, xi, Xi); end % 第二个非零项右侧基函数贡献 w2 0; den2 Xi(ip1) - Xi(i1); if den2 0 w2 (Xi(ip1) - xi) / den2 * bspline_basis(i1, p-1, xi, Xi); end N w1 w2; end基函数计算是整个程序的地基对每个高斯积分点都要调用。注意代码里两个分母判断了大于0才计算这是必须的节点向量里允许出现重复节点重复时对应区间长度为0如果直接除会得到NaN一传十十传百最后整个刚度矩阵报废。传入的Xi必须是完整节点向量两端节点要重复p1次这是让边界基函数具备插值性的前提。如果你想验证代码对不对可以在MATLAB里画一下p2、Xi[0 0 0 1 1 1]时的三条基函数曲线它们应该是对称的山丘形。2.2 权重与有理化NURBS基函数的完整计算只用B样条的话圆、圆柱面这类曲线曲面表达不出来等几何分析最值钱的能力就废了一半。NURBS在B样条基础上给每个控制点配一个权重基函数变成有理形式。权重参与计算的代码可以这样组织function [R, dRdxi] nurbs_basis(i, p, xi, Xi, w) % 计算带权重的NURBS基函数值及其一阶导数 % w: 所有权重组成的向量与基函数数量一致 % 先算无权重B样条基函数值 N bspline_basis(i, p, xi, Xi); % 算分母所有基函数加权求和 Wsum 0; for k 1:length(w) Wsum Wsum w(k) * bspline_basis(k, p, xi, Xi); end % 有理基函数值 R w(i) * N / Wsum; % 导数用商的求导法则 dN bspline_basis_deriv(i, p, xi, Xi); dWsum 0; for k 1:length(w) dWsum dWsum w(k) * bspline_basis_deriv(k, p, xi, Xi); end dRdxi (w(i) * dN * Wsum - w(i) * N * dWsum) / Wsum^2; end这里的bspline_basis_deriv是B样条基函数的一阶导由递推关系写成function dN bspline_basis_deriv(i, p, xi, Xi) if p 0 dN 0; return; end term1 bspline_basis(i, p-1, xi, Xi) / (Xi(ip) - Xi(i)); term2 bspline_basis(i1, p-1, xi, Xi) / (Xi(ip1) - Xi(i1)); dN p * (term1 - term2); end导数计算里我用den0判断而不是deneps因为节点重复时分子也同时为00/0这个未定式在MATLAB里产生NaN但0/eps也还是0乘0不影响结果。nurbs_basis这个函数的输出R和dRdxi是后续雅可比矩阵计算的直接输入权重Wsum在几何映射里还要再除一次这点容易漏如果几何映射里已经除了权重形函数计算里就不该再除一遍否则权重被平方了几何完全变形。2.3 张量积基函数与几何映射从参数域到物理域二维问题里NURBS基函数是张量积结构二维基函数等于两个一维基函数的乘积。这套程序里几何映射函数按这个结构组织function [x, y] nurbs_geometry_point(xi, eta, P, W, Xi, Eta, p, q) % 将参数域点(xi, eta)映射到物理坐标(x, y) % P: (n1)x(m1)x2 控制点坐标矩阵 % W: (n1)x(m1) 控制点权重矩阵 % p, q: 两个参数方向的基函数阶数 x 0; y 0; denom 0; for i 1:size(P,1) for j 1:size(P,2) % 张量积基函数乘以对应权重 R nurbs_basis(i, p, xi, Xi, W(i,:)) ... * nurbs_basis(j, q, eta, Eta, W(:,j)) / Wsum_total; denom denom R; x x R * P(i,j,1); y y R * P(i,j,2); end end x x / denom; y y / denom; end注意这个循环里的Wsum_total要单独算不能直接复用前面nurbs_basis里的局部分母。正确做法是先把所有非有理B样条基函数乘以权重求和得到全局分母再统一归一化。一个容易踩的坑是如果每个方向调用nurbs_basis时都除了一次该方向的权重和再在张量积里除一次总权重权重会被重复归一化。我在第一次调试这套程序时就在这里翻过车最后用单位圆几何做测试如果映射出的点全部偏离半径1多半就是权重除多了。几何映射的输入输出关系是参数域是一个单位矩形物理空间是这个矩形经过控制点和权重扭曲后的形状。程序里所谓“网格”其实只是控制点连成的控制多边形真正的几何是所有高斯积分点映射后的坐标。理解这一层后面看刚度矩阵就顺了。3. 单元刚度矩阵与全局组装从高斯积分到系统方程3.1 单元刚度矩阵雅可比矩阵与B矩阵的构建刚度矩阵的数学形式和传统有限元完全一致差别只在形函数和导数来源不同。二维线弹性平面应力问题的单元刚度矩阵表达式是Ke ∫ B^T D B |J| dξ dη其中D是弹性矩阵由弹性模量和泊松比决定。这套程序在单元刚度计算时做的是三层循环先遍历高斯点再在积分点处组装B矩阵最后累加。核心片段如下function Ke element_stiffness(elem, controlPts, weights, Xi, Eta, p, q, D) % elem: [i j] 节点向量中的单元索引 % controlPts: 该单元关联的局部控制点坐标 % 返回: 单元刚度矩阵 [gp_xi, gw_xi] gauss_rule(p1); [gp_eta, gw_eta] gauss_rule(q1); n_local (p1) * (q1); Ke zeros(2 * n_local); for k 1:length(gp_xi) % 把高斯点从[-1,1]映射到节点区间 xi (Xi(elem(1)) Xi(elem(1)1)) / 2 ... gp_xi(k) * (Xi(elem(1)1) - Xi(elem(1))) / 2; jac_xi (Xi(elem(1)1) - Xi(elem(1))) / 2; for l 1:length(gp_eta) eta (Eta(elem(2)) Eta(elem(2)1)) / 2 ... gp_eta(l) * (Eta(elem(2)1) - Eta(elem(2))) / 2; jac_eta (Eta(elem(2)1) - Eta(elem(2))) / 2; % 局部基函数值和参数域导数 [R, dR_dxi, dR_deta] local_basis_2d(elem, controlPts, weights, xi, eta); % 雅可比矩阵参数坐标到物理坐标的转换 J [sum(dR_dxi .* controlPts(:,1)), sum(dR_dxi .* controlPts(:,2)); sum(dR_deta .* controlPts(:,1)), sum(dR_deta .* controlPts(:,2))]; detJ det(J); invJ inv(J); % 物理空间导数 dR_dx invJ(1,1) * dR_dxi invJ(1,2) * dR_deta; dR_dy invJ(2,1) * dR_dxi invJ(2,2) * dR_deta; % 二维平面应力B矩阵每个控制点两个自由度 B zeros(3, 2*n_local); B(1, 1:2:end) dR_dx; B(2, 2:2:end) dR_dy; B(3, 1:2:end) dR_dy; B(3, 2:2:end) dR_dx; w gw_xi(k) * gw_eta(l) * jac_xi * jac_eta * detJ; Ke Ke B * D * B * w; end end end这里最值得说明的是B矩阵的排列方式每个控制点有x、y两个自由度所以B矩阵的列数是局部控制点数的两倍奇数位放x方向的导数、偶数位放y方向的导数。这样排的目的是一旦组装完Ke直接把它塞进全局矩阵自由度顺序与位移向量U的排列一致U [u1_x; u1_y; u2_x; u2_y; ...]不需要再调整列顺序。雅可比矩阵的行列式detJ是参数域微小面积到物理域微小面积的缩放因子必须乘进积分权重它的值恒正说明控制点排列方向正确一旦出现负值说明控制点索引顺序反转了需要换方向。3.2 高斯积分点数选择精度与组装时间的权衡这套程序默认每个方向取p1个高斯点也就是二次基函数每方向取3个高斯点。这个数量对等几何分析来说是被实践反复验证过的经验值既能保证刚度矩阵在多项式部分完全精确积分又不至于把组装时间拖长。需要说明的是NURBS引入有理分式后严格来说p1个高斯点不能保证被积函数完全精确因为有理分母破坏了多项式性质。但实际算下来误差在10的负7次方量级完全够工程用。如果你想提高积分精度盲目把高斯点从3加到6效果不明显但组装时间几乎翻倍。我之前做过一次对照同样一个4000自由度的曲梁算例p2时n3和n6的位移解在小数点后第七位才出现差异组装时间增加了45%。自由度没变精度也没实质变化这钱花得冤枉。正确的调参思路是先保证网格单元数量足够再考虑加高斯点。等几何分析的收敛性主要由基函数阶数和节点加密策略决定增加高斯点只是消除积分误差并不能替代h加密和p加密。3.3 全局组装局部控制点到全局自由度的索引映射单元刚度矩阵算出来之后麻烦事才真正开始。全局刚度矩阵的尺寸是2×控制点总数每个单元关联的控制点不像传统有限元那样恰好是单元顶点或者边中点而是由节点向量跨度决定的局部控制点块。把局部索引映射成全局索引是组装环节最容易写错的地方。% 组装全局稀疏矩阵 n_dofs 2 * size(controlPts, 1); K sparse([], [], [], n_dofs, n_dofs, 20 * n_dofs); F zeros(n_dofs, 1); for e 1:size(active_elems, 1) Ke element_stiffness(active_elems(e,:), controlPts, weights, Xi, Eta, p, q); % 局部自由度到全局自由度的编号 local_dofs element_dofs(active_elems(e,:), p, q); K(local_dofs, local_dofs) K(local_dofs, local_dofs) Ke; end用稀疏矩阵存储是性能关键。如果你第一版代码用全稠密矩阵写几千自由度就顶不住了而用sparse预分配内存几万自由度也就几秒的事。另一个技巧是不要用sparse(zeros(n,n))这种写法它会先在内存里建一个巨大的稠密零矩阵再转稀疏浪费时间和内存。直接按sparse(i,j,v,m,n)用三元组方式构造或者用spalloc预分配非零元素数量速度会快一个数量级。element_dofs函数的实现逻辑是先算出当前单元在节点向量中的起始索引再根据p和q阶数推导出关联的控制点行和列范围。这个映射关系一旦错位组装出来矩阵的带宽分布就会混乱常见表现是刚度矩阵秩亏或者位移场一片混乱。4. 边界条件与求解从方程组到可视化结果4.1 强制位移边界条件行替换法的实现与取舍等几何分析里的边界条件施加比传统有限元要多一步思考。NURBS基函数在边界上的控制点确实对应几何边界但物理意义和插值节点不同边界控制点的位移值不是直接等于边界位移而是通过边界基函数加权后得到边界位移。好在这套程序处理的是二维简单模型边界控制点沿边界方向只有一重基函数可以直接对控制点自由度施加位移。常用的是行替换法% 施加强制边界条件 for i 1:length(boundary_dofs) dof boundary_dofs(i); val prescribed_disp(i); % 将已知自由度移动到方程组右侧 F F - K(:, dof) * val; % 清零对应行列后对角线置1 K(dof, :) 0; K(:, dof) 0; K(dof, dof) 1; F(dof) val; end这段代码的好处是它对刚度矩阵的对称性没有破坏因为行列同时清空留下的系统仍然是正定的可以直接用K\F求解。它的代价是如果你用迭代法解方程这种修改方式会破坏原来的预处理信息但在自由度规模几千到一万这个区间直接求解器完全不是问题。另一种常见做法是配点法在边界上选取物理点要求这些点处的位移精确等于给定值再把约束方程叠加到系统矩阵。这种方法的优势是不需要修改刚度矩阵的对角元素但对边界点的选择比较敏感选少了约束不足选多了方程过约束。对一般的平面问题行替换法简单直接我建议你第一版先用它跑通再去研究配点法。4.2 全局求解与自由度重排从一维解向量到物理位移场全局方程组装好之后求解就是一行代码的事U K \ F;MATLAB会自动识别K是稀疏对称正定矩阵走LDL分解路径。如果模型自由度超过五万可以考虑换用ichol预处理共轭梯度法但普通学习和论文复现场景下K\F完全够用。求解之后得到的U是一个长度为2×控制点数的列向量需要重组成分量的形式才能做后处理ncp size(controlPts, 1); u_field zeros(ncp, 2); u_field(:,1) U(1:2:end, 1); u_field(:,2) U(2:2:end, 1);控制点位移场不等于物理点的最终位移要得到任意物理坐标处的精确位移还需要通过几何映射插值。这一步在等几何分析里是顺理成章的既然位移场用同一套NURBS基函数表达那对任何参数坐标点(xi, eta)位移就等于局部控制点位移乘以局部基函数再求和。这个性质让IGA的后处理精度远高于传统的“节点位移等参插值”。4.3 可视化与误差评估控制网格、变形场与应变云图这套程序的后处理模块能直接绘制三样东西控制网格、变形后几何、应变云图。绘制控制网格的方法是连接控制点能直观看到基函数支撑和几何形态绘制变形几何则需要把控制点位移加到控制点坐标上重新做几何映射。% 在网格线上密集采样绘制变形后的几何 nPlot 40; [xg, yg] meshgrid(linspace(0,1,nPlot), linspace(0,1,nPlot)); xd zeros(size(xg)); yd zeros(size(yg)); for r 1:nPlot for s 1:nPlot [x_orig, y_orig] nurbs_geometry_point(xg(r,s), yg(r,s), controlPts, weights); [dx, dy] nurbs_displacement_point(xg(r,s), yg(r,s), u_field, weights); xd(r,s) x_orig scale * dx; yd(r,s) y_orig scale * dy; end end surf(xd, yd, zeros(nPlot), DisplayName, Deformed Geometry);参数说明scale是位移放大倍数因为实际位移通常只有几何尺寸的千分之一不放大根本看不出变形形态。nPlot是采样密度越大曲线越光滑但计算时间线性增长通常取40到60就已经足够。判断程序算得对不对的简单办法是选一个已知解析解的模型比如悬臂梁端部受集中力把端部位移和理论值对比。误差在1%以内说明基函数、积分、边界条件基本没有大问题。5. 排查与避坑IGA程序最常见的四个失败点5.1 递归除零导致整套矩阵变成NaN现象是位移解全部输出NaN甚至连第一步组装时Ke就出现NaN。原因是节点向量里有重复节点B样条基函数递归时某一段定义区间长度为0没有保护直接触发了0除以0。排查方法是在基函数函数入口加一个assert(isfinite(N))第一个出现NaN的基函数序号会立刻暴露问题节点。解决办法是在所有分母处加上长度大于eps的判断条件分子为0时跳过该分支。我从这套程序里学到的习惯是只要写了除法分母一律考虑“可能是0”的情况不要依赖数据天然正常。5.2 形函数导数做数值差分应力场出现锯齿现象是位移解正常但应力场分布出现规律性锯齿误差集中在单元边界附近。原因是某个版本里为了省事对形函数导数用了有限差分近似dN (N(xih)-N(xi-h))/(2h)。这种做法在高斯点处计算导数的精度被截断误差限制而刚度矩阵的B矩阵对导数误差非常敏感最后污染的是整个应力云图。解决方法是恢复解析递推公式数值差分只能用于校验正确性绝不能用进组装流程。从那以后我每次计算导数都会做一个一致性检查把数值差分算出来的导数和解析递推结果对比误差应该在机器精度级别如果出现10的负6次方量级的差异一定是递推公式某个分支写错了。5.3 雅可比行列式出现负值变形方向反转现象是应变能出现负值或者施加向下载荷结果模型向上变形。原因是控制点矩阵的排列方向与参数域方向不匹配导致雅可比行列式在整个单元域上恒为负值。等几何分析中的雅可比行列式描述的是从参数域到物理域的面积映射它的正负决定了应力和应变的正负约定是否一致。排查方法是在单元刚度函数里加一行打印最小detJ如果小于某个阈值比如1e-10就说明控制点顺序有问题。解决方法是把控制点数组按x坐标升序、y坐标升序重新排列同时调整节点向量索引关系。这套程序里我吃过一次亏建立控制点矩阵时用了顺时针排序结果整个模型的应力场全部反号怎么看怎么不对劲。5.4 局部-全局索引错位刚度矩阵奇异现象是求解阶段报Matrix is singular或者位移场出现完全不符合物理直觉的剧烈突变。原因是组装单元刚度矩阵时局部控制点编号到全局自由度的映射错了一位。NURBS单元的关联控制点不是按传统单元的节点顺序而是由节点向量跨度决定的一个子矩阵块写映射时经常出现差一错误。排查方法是组装完成后随机抽某个单元把它局部自由度和全局自由度编号打印出来人工核对前两行非零元素的位置是否落在预期的列区间内。更稳妥的办法是写单元级测试对一个单单元模型手动推导二乘二控制点的刚度矩阵前几项和程序输出逐项对比全对再进下一阶段。索引类问题不会在某个算例里突然出现它一旦出现就是系统性的所以值得花半天时间把映射函数验证到滴水不漏。6. 二次开发思路从二维线弹性走向真正的实用6.1 几何非线性扩展替换本构和构型即可起跑原程序做的是小变形线弹性改成大变形要动的核心不是基函数而是“构型更新”刚度矩阵的B矩阵要换成当前构型下的形式同时引入格林-拉格朗日应变和柯西应力。具体做法是在每次迭代中把控制点坐标加上当前位移增量重新计算雅可比和B矩阵再组装切线刚度矩阵。原有基函数、几何映射和组装框架都可以原封不动复用主要工作量在写非线性求解器那一层。这个改造路径充分体现了等几何分析的优势几何每步都在变但基函数和参数域不变所有中间量都在同一套NURBS框架里计算。6.2 多面片耦合与拓扑优化从单位圆到复杂装配体原程序假设整个模型是单张NURBS面片实际工程几何几乎都是多面片拼起来的。跨面片模型要把交界面处控制点自由度耦合起来常见做法是加罚函数约束或者拉格朗日乘子法。改造时先在程序中加一个“面片交界面检测”模块输出重叠面上的控制点索引对再在全局组装时把约束项拼进刚度矩阵。另一个高价值方向是拓扑优化以控制点坐标为设计变量把灵敏度和基函数结合就能在CAD几何上直接做密度场拓扑优化绕开传统优化里反复网格重构的循环。这就是等几何分析在当前学术和工业里最活跃的两条应用线。6.3 保持代码可扩展性的一点习惯我想给拿到这套程序准备二次开发的人一个建议保持基函数计算、几何映射、单元刚度三个模块的接口纯函数化——输入参数坐标输出值和导数绝不读写全局变量。我第一次做扩展时图省事在刚度矩阵里直接引用了全局控制点数组结果后续加多面片耦合时被迫把所有函数签名重写了一遍耗时整整一周。抛开这个故事不谈纯函数接口还能让你在单元级别做独立测试每条曲线、每张曲面都能先验证后装配。从那以后我每次做IGA相关开发都强制走一遍“先单元级验证再全局组装”的流程基函数测一遍、几何映射测一遍、刚度矩阵测一遍全部见过绿灯才继续扩展。希望帮到你。本文还有配套的精品资源点击获取