简介：面向结构稳定分析的 MATLAB 弧长法实现脚本，适合需要处理非线性屈曲路径与临界荷载计算的结构工程师、研究人员和高年级学生。压缩包内共 2 个 m 文件（Arclength.m 与 Arclength2.m），整体大小约 5KB，分别对应弧长法的基本求解框架与扩展分析版本，可帮助使用者快速搭建屈曲分析流程。已有 500 人学习下载。脚本围绕虚拟弧长参数控制步长、非线性方程组迭代、几何与材料非线性等核心内容展开，既有基础实现也有改良思路，便于对照学习屈曲荷载识别、后屈曲路径追踪以及多自由度系统扩展。通过理解这两个脚本，读者能掌握弧长法在 MATLAB 中的编程要领，并直接应用于结构稳定性的初步分析或二次开发。
1. 弧长法不是黑匣子：结构稳定分析中那条追不回的后屈曲路径
做结构稳定分析的人，多半在某个晚上遇到过同一件事：一条挺漂亮的荷载-位移曲线，算到极值点附近，牛顿迭代突然发散，屏幕上一片 NaN。你试着把荷载增量调小，曲线还是断在峰值前。这时候基本可以确定，问题不在网格、不在材料，而在你用的是荷载控制还是位移控制，且这两者在后屈曲路径跟踪里都有硬伤。弧长法（arc-length method）就是用来解决这个问题的：把荷载因子也当成未知数，让每一步迭代沿着一条弧长约束前进，跨过极值点，继续描出完整的后屈曲路径。这个思路在 b幈ckling 分析里几乎是标配，尤其薄壁结构、拱结构、壳结构的稳定分析。
标题里的 Arc-length.rar 这类资源包，网上常年能搜到，里面核心脚本无非是 Newton-Raphson 的扩展版加一两个算例。真正难的不是几十行迭代代码，而是预测步方向、λ 符号维护、选根策略、弧长自适应这几个细节。这篇笔记把这些拆开讲清楚，给出可以直接改到自己程序里的 MATLAB 实现。适合正在做有限元二次开发、写研究生论文、或者被后屈曲路径折磨的工程师。
2. 为什么非线性稳定分析离不开弧长法：极值点附近的牛顿法为何翻车
2.1 荷载控制与位移控制在失稳分析里的局限
常规非线性静力分析，最常见的是荷载控制下的 Newton-Raphson 迭代。每个子步外荷载增量固定，迭代过程中只更新位移。这个做法在一般弹塑性问题里很稳，但一碰极限点就失效。原因很简单：在极限点处，结构切线刚度矩阵奇异，荷载增加不再对应唯一的位移增量，甚至结构需要卸载才能维持平衡。你给定了一个正的荷载增量，系统却要求负的荷载增量，迭代自然发散。这不是增量步长取大了，而是控制方程本身就无解。
位移控制能部分解决极值点问题。把加载方式改成指定位移增量，相当于换了一个控制参数，切线刚度矩阵奇异问题被绕开。但位移控制只对“荷载随位移单调变化”的情况有效。一旦遇到回跳型失稳，也就是 snap-back，同一个位移可能对应三个荷载状态，位移控制也会失去唯一的投影方向。工程结构里这类情况并不罕见，浅拱、扁壳、负高斯曲率曲面都可能出现。
这两种控制方式共同的特点是：控制变量是单一标量，要么是荷载，要么是位移。而弧长法的出发点是放弃这种单一控制，把荷载因子 λ 和位移增量一起放进未知数，用一条弧长约束把两者绑定。这样一来，极值点和回跳点都只是路径上的普通点，不再需要特殊处理。
2.2 三种弧长形式：柱面、球面、线性弧长怎么选
弧长法按约束方程的不同分成几种。最经典的是 Riks 提出的线性弧长法，约束方程里位移增量和荷载增量成线性关系。后来 Crisfield 做了改进，提出球面弧长法，约束方程中同时含位移增量和荷载增量的二次项。工程实现里最常用的是柱面弧长法，它只约束位移增量向量的范数，荷载因子不再出现在约束方程中。
柱面弧长法之所以用得多，是因为它在实现上最省心：每步迭代只需要维持位移增量的长度不变，二次方程中的二次项只和位移相关，数值稳定性比球面弧长法好控制。球面弧长法理论上对 snap-back 路径更鲁棒，但二次方程中可能出现 a 接近零的情况，程序处理起来比较麻烦。线性弧长法现在用得越来越少，因为它的切平面近似在后屈曲路径上误差偏大，往往需要更小的步长才能跟踪。
选型建议很简单：先写柱面弧长法，绝大多数结构稳定分析问题都够用。真遇到球面弧长才能收敛的算例，再改约束方程也不迟。两者的主循环几乎一样，差别只在二次方程的系数上。
2.3 弧长法的两个方程与迭代路径
弧长法在每个增量步内要同时满足两个方程。第一个是平衡方程：
r(u, λ) = f_int(u) - λ f_ref = 0
第二个是柱面弧长约束方程：
Δu^T Δu = Δl^2
f_int 是内力向量，f_ref 是参考荷载向量，λ 是荷载因子，Δu 是当前子步从起点开始的累计位移增量，Δl 是当前弧长半径。每一步迭代的目标，是找到一组 (u + Δu, λ + Δλ)，让这两条方程同时成立。
整个迭代过程分预测和校正两步。预测步从当前切线刚度出发，沿着切线方向走一个弧长半径，得到一个初始猜测点。校正步在这个猜测点附近做 Newton 迭代，每次求解都需要用到块消元，因为荷载因子是额外引入的未知数，不能像常规非线性分析那样只求解位移增量。迭代收敛后，检查当前子步是否满足平衡和弧长约束，不满足就继续校正，满足就进入下一个子步，并重新计算下一弧长半径。
这套流程听起来不复杂，真正写出代码也就几十行。但每个环节都有细节，下面一章直接给出实现。
3. 用 MATLAB 手写一个弧长法求解器：预测步、校正步与块消元
3.1 预测步：切线位移与 λ 增量的符号维护
每个增量步开始时，先用当前切线刚度矩阵 K 求参考荷载产生的切线位移：
% 子步开始：用当前切线刚度做预测 dut = K \ f_ref; % 参考荷载下的切线位移 lam_inc = sign_l * dl / sqrt(dut' * dut); % 荷载因子增量 du = lam_inc * dut; % 预测位移增量 lam = lam + lam_inc; u = u + du;这里的 K 是当前状态的切线刚度矩阵，f_ref 是参考荷载向量。dut 的物理含义是单位参考荷载下结构会往哪个方向走，它的范数大小反映了结构在当前状态的柔度。lam_inc 的绝对值由弧长半径 dl 除以 dut 范数得到，因为预测点的位移增量长度要等于 dl。sign_l 是符号标志量，记录上一子步 λ 增量的符号。
注意，预测步不迭代，它只负责给出一个可靠初值。如果预测方向错了，后面校正步再努力也可能翻车。符号维护的第一版实现，用一个全局变量记录上一子步的 λ 增量符号，第一个子步根据参考荷载方向给正号。后屈曲路径上这个符号可能翻转，如果翻转发生在极值点前，预测就会反向。后面避坑章会讲更稳的方向判断方式，这里先按最简单的符号记忆处理。
3.2 校正步：块消元求 λ 增量与选根
校正步的核心，是把每个 Newton 迭代步的位移修正拆成两部分：一部分由不平衡力引起，另一部分由荷载因子增量引起。然后用弧长约束方程解出 λ 增量。这段是弧长法最容易写错的地方。
for iter = 1:max_iter [K, r] = assemble(u, lam); % 组装切线刚度与平衡残差 du_r = K \ (-r); % 不平衡力引起的位移修正 du_t = K \ f_ref; % 单位荷载因子引起的位移修正 a = du_t' * du_t; b = 2 * du_t' * (du + du_r); c = (du + du_r)' * (du + du_r) - dl^2; disc = b^2 - 4 * a * c; if disc < 0 dl = dl * 0.5; % 半径过大，减半后重来 u = u_start; lam = lam_start; du = du_start; continue; end sq = sqrt(disc); dlam1 = (-b + sq) / (2 * a); dlam2 = (-b - sq) / (2 * a); % 选根：取使增量方向与上一步增量方向内积更大的根 cos1 = (du + du_r + dlam1 * du_t)' * du; cos2 = (du + du_r + dlam2 * du_t)' * du; if cos1 > cos2 dlam = dlam1; else dlam = dlam2; end du = du + du_r + dlam * du_t; u = u + du_r + dlam * du_t; lam = lam + dlam; if norm(r) < tol * norm(f_ref) % 力残差收敛判据 break; end end这段代码里的 assemble 函数是占位符，实际项目中替换成你自己的单元刚度组装和内力计算函数。a、b、c 是二次方程的系数，这个二次方程来自弧长约束：把 du + du_r + dlam * du_t 代进 Δu^T Δu = dl^2，展开后就是 a * dlam^2 + b * dlam + c = 0。两个根都满足弧长约束，但只有一个根对应真实的平衡路径，选根要从几何上判断：计算两种候选增量方向与当前子步已有增量方向的内积，取内积更大、也就是方向更接近的那个根。
判别式 disc 小于零，说明当前弧长半径下，约束圆和平衡路径不相交。出现这种情况，最常见的处理是减半半径之后回到子步起点重新算，而不是在当前点上硬凑。u_start、lam_start、du_start 要在子步进入时保存，作为重试的后悔药。
3.3 能跑的验证算例：两杆桁架的 snap-through
弧长法写完后，不能直接拿去算复杂模型，先用一个教科书级算例验证程序逻辑。两杆桁架的 snap-through 问题是首选：结构简单，只有两个自由度，且对称约束后只有一个竖向位移，但它的荷载-位移曲线具有完整的极值点、下降段和二次上升段，能检验弧长法的每一步。
function [K, r] = truss_state(y, EA, L, h, lam, P) l0 = sqrt(L^2 + h^2); % 初始杆长 l = sqrt(L^2 + (h - y)^2); % 变形后杆长 N = EA * (l - l0) / l0; % 轴力，压缩为负 f_int = 2 * N * (h - y) / l; % 杆件给节点的合力（向上为正） r = -f_int - lam * P; % 平衡残差，外载向下为正 K = 2 * EA / l0 * (1 - l0 * L^2 / l^3); % 切线刚度 end这个函数返回的 K 是标量，因为对称性让水平位移始终为零。外荷载 P 取 1000 N，EA 取 1e5 N，L 取 0.5 m，h 取 0.1 m。把 assemble 函数换成 truss_state，主循环不变，就能得到一条完整的 snap-through 曲线：荷载先随位移上升，到达第一个极值点后下降，结构跳到下稳定分支，然后继续上升。这就是弧长法的标志性能力，它让牛顿迭代在下降段也能收敛。
验证时留意一点：如果程序在极值点附近依然发散，先不要怀疑弧长法，检查选根逻辑。打印出每一步的两个内积值，通常能看出选根选反了。这个算例跑通后，再把它替换成你自己的梁单元、壳单元或者实体单元，弧长法主循环一行都不用改。
3.4 把弧长法嵌进你自己的有限元程序
我自己在项目里的组织方式，是维护三个独立模块。第一个是状态函数，输入节点位移和荷载因子，输出切线刚度矩阵和残差向量，这是和你单元库唯一相关的部分。第二个是弧长法主循环，只调用状态函数，不关心单元类型。第三个是后处理脚本，负责提取荷载-位移曲线、更新弧长半径、输出增量步中间结果。
模块划分决定了调试效率。新手最容易犯的错是把弧长法逻辑和单元组装修在一起，最后程序跑不起来时，分不清是几何非线性有问题，还是选根有问题。我一般让状态函数先和普通牛顿法配合，确认单点加载的弹塑性分析能收敛后，再套上弧长法外壳。这样一旦后屈曲路径出错，问题基本锁定在弧长相关代码里。
4. 弧长法的参数怎么给：初始半径、自适应与收敛容差
4.1 初始弧长半径的经验取值
初始弧长半径 dl 是弧长法里影响最大的参数。给太大，第一步就会越过极值点，二次方程判别式为负，程序不断减半重试，效率极低。给太小，整个计算步数太多，一个复杂模型可能要跑几百步。常见做法是先跑半步切线预测，看参考荷载会产生多大的位移：
% 估算初始弧长半径 du_ref = K0 \ f_ref; dl0 = 0.1 * sqrt(du_ref' * du_ref);0.1 这个系数是我常用的起点。它表示第一增量步的位移长度大概是参考荷载静力位移的十分之一。结构较软、参考荷载取得偏大时，需要把这个系数调小到 0.01；结构很硬、路径简单，可以给到 0.2 甚至 0.5。注意 dl0 的量纲是位移的范数，和你模型的单位制直接相关。如果用 mm 建模，dl0 就是多少毫米；如果用 m 建模，就是多少米。
另一个经验是先用线性屈曲分析估算临界荷载，再把参考荷载取到临界荷载的 1.2~2 倍。这样 dl0 的数值处于一个合理的位移量级，不至于出现“参考荷载太大，切线位移范数跑到几百毫米”的情况。
4.2 自适应弧长与收敛容差
固定弧长半径能跑，但效率不高。后屈曲路径上曲率变化剧烈，固定半径会导致在极值点附近反复减半，在平坦段又浪费步数。自适应弧长的标准做法是按上一子步的迭代次数调整半径：
% 自适应弧长：目标迭代次数 n_target，实际迭代次数 n_iter dl = dl * sqrt(n_target / n_iter); dl = min(max(dl, dl_min), dl_max);这个公式的逻辑是：如果某个子步只用了两三步就收敛，说明路径比较平缓，下一步可以走得更远；如果用了十步才收敛，说明曲率大，下一步要收紧。目标迭代次数我一般设 5，下限 dl_min 取 dl0 的 0.05 倍，上限 dl_max 取 dl0 的 5 倍。限幅必须加，否则平坦段会把半径放大到离谱，遇到斜率突变时又来不及收回来。
收敛容差方面，不要只用位移增量判据。极值点附近，位移对荷载的变化率趋于无穷，位移判据会给出"已经收敛"的错误信号。我一般用力残差判据，norm(r) 小于 tol 乘以 norm(f_ref)，tol 取 1e-6。荷载因子也参与迭代了，理论上应该检查 λ 增量是否进入容差，但实际工程中力残差结合最大迭代次数限制已经足够稳定。
4.3 配合参考荷载与材料非线性时的参数注意
参考荷载向量 f_ref 的选择比很多人想象的更关键。它不一定要等于真实荷载，只是一个方向向量，真实荷载由 λ f_ref 给出。但如果 f_ref 的分布形态和真实荷载差距过大，比如真实荷载是集中力、f_ref 却给成均布力，后处理的荷载-位移曲线会很难看懂。我一般按真实荷载的分布形态填 f_ref，大小取预估极限荷载。
考虑材料非线性时，弧长法主循环不用改，但状态函数里要用材料切线刚度。注意一点：材料进入塑性后，卸载路径和加载路径的模量不同，弧长法在跨越极值点时可能会进入卸载，这时如果状态函数里没区分加卸载，结果曲线会失真。这个问题在理想弹塑性模型里特别明显。
几何非线性是大前提。弧长法处理的是几何失稳，状态函数里的应变-位移关系必须包含大变形项。如果单元还是小变形假设，弧长法参与迭代的切线刚度矩阵是常数，后屈曲路径根本不存在。
5. 弧长法避坑指南：四个高频翻车现场与排查思路
5.1 增量步预报方向错误导致“往回跑”
现象：荷载-位移曲线在极值点前就开始往回走，或者曲线整体沿加载反方向展开，后处理里看到的是一条镜像路径。
原因：预测步的符号标志 sign_l 沿用上一子步，但路径在某个极点发生方向翻转，符号没有跟着变。这通常发生在极值点附近、弧长半径偏大时。
解决：不要把 sign_l 只存上一个子步的符号。我常用的做法是在每个子步预测前，计算当前切线位移 du_t 和上一子步总位移增量 du_old 的内积，如果内积为负，说明路径方向要翻转，强制 sign_l 取反。这个判断在大多数工程问题里足够可靠。
5.2 切线刚度阵奇异导致线性求解漂移
现象：主循环迭代过程中，K 的条件数爆掉，MATLAB 给出警告，du_r 或 du_t 出现巨大数值，曲线直接飞出去。
原因：到达理论极限点时切线刚度矩阵奇异，这是必然的。MATLAB 的 A\b 在矩阵接近奇异时会给出带警告的数值解，但误差已经被放大，校正步不可能收敛。
解决：一是诊断，把 cond(K) 打印出来，看奇异发生的位置是否和极值点重合。二是处理，在 K 上叠加一个小扰动：K + 1e-8 * norm(K) * eye(ndof)。这个技巧不改变路径形态，只是让线性求解器不炸。三是策略调整，遇到奇异时强制减半弧长半径，让每个子步的增量变小，避开精确的奇异点。记住，奇异是结构性质，不是程序 bug，弧长法只是尽量跨过它，不能在奇异点上求逆。
5.3 二次方程判别式一直为负，弧长减半死循环
现象：程序进入减半重试流程后，半径一减再减，直到低于 dl_min 还是 disc < 0，最终停在某一步不再前进。
原因：除了半径过大，另一个常见原因是预测方向垂直于真实路径。柱面弧长法要求预测点落在约束圆附近，如果预测方向完全偏掉，约束圆和路径就是不交，减半半径只能让情况慢慢好一点，但可能永远到不了可接受范围。
解决：判断不仅是减半半径，还要重置预测方向。我一般会强制从当前的 K 重新生成 du_t，然后让符号标志翻转一次，给一个反向预测。这种“试探性反向”在分叉点附近特别有效。同时检查 K 是否已经奇异，如果 cond(K) 超过 1e12，先把 K 扰动稳住，再做反向预测。
5.4 结果比实验值高一大截：初始缺陷被丢掉了
现象：弧长法算出来的失稳荷载是实验值的好几倍，曲线形态也不对，实验里明显的下降段算出来却是一条平缓上升线。
原因：理想结构的屈曲是分支点型失稳，弧长法追踪的平衡路径从完美几何出发，而真实结构永远有初始几何缺陷、偏心荷载和残余应力。忽略缺陷等于求解另一个结构。
解决：这是弧长法本身的问题，不能靠调参数解决。正确流程是先做线性屈曲分析，拿到一阶屈曲模态，再把模态乘以一个比例系数叠加到初始几何上，然后重新跑弧长法。这个流程几乎是所有稳定分析的标准开头。具体做法放在下一章，属于把弧长法用透的进阶操作。
6. 把弧长法用透：特征值屈曲预判临界荷载与缺陷敏感性
6.1 先算线性 buckling，拿到临界荷载和模态
直接上弧长法之前，花几分钟算一次线性屈曲分析，收益很高。线性屈曲的特征值问题写出来是：
(K0 + λ Kg) v = 0
K0 是小变形切线刚度，Kg 是几何刚度矩阵，λ 是屈曲因子。在 MATLAB 里用 eig 求解广义特征值问题：
% 线性屈曲：求解广义特征值问题 [V, D] = eig(K0, -Kg); [eigen, idx] = sort(diag(D)); lambda_cr = eigen(1); % 最小特征值为临界荷载因子 phi = V(:, idx(1)); % 对应一阶屈曲模态这里把 -Kg 作为第二参数，特征值直接就是屈曲荷载因子。取 eigen 的最小值，对应的特征向量就是结构最弱的失稳形态。这一步的意义不只是给弧长法一个参考荷载，更重要的是为下一步构造初始缺陷提供模态。
6.2 用一阶屈曲模态构造初始几何缺陷
拿到模态 phi 后，把它归一化，再乘一个缺陷幅值 alpha，叠加到节点坐标上：
% 叠加初始几何缺陷：alpha 取跨度的 1/1000~1/300 u_init = u_init + alpha * phi / max(abs(phi));alpha 的取值按工程惯例来。钢结构一般取构件跨度的 1/1000，薄壳结构偏保守取 1/300。如果你手头有施工规范或实验数据，按规范给缺陷幅值，没有的话就从 1/1000 起步试。
这个带缺陷的模型再跑弧长法，得到的荷载-位移曲线才是工程上真正关心的后屈曲响应。缺陷幅值对极限荷载的影响程度，就是结构的缺陷敏感性。改成不同 alpha 跑一遍，把极值点荷载画成一条曲线，就能看出这个结构是缺陷敏感型还是不敏感型。这个过程几乎不用改弧长法代码，只是换初始几何。
6.3 从曲线读极限承载力与失稳类型
弧长法输出的是完整的荷载-位移曲线。带缺陷模型的曲线第一个极值点对应的荷载因子乘以参考荷载，就是考虑初始缺陷后的极限承载力。如果曲线在峰值附近突然下降，说明是极值点失稳；如果曲线在某个荷载下从主线跳向支线，则可能伴随分叉点失稳。
我现在的习惯是：任何结构稳定分析，先线性 buckling 估算临界荷载，再按规范比例叠加一阶模态，最后用弧长法复算后屈曲路径。这个流程跑通了，无论多复杂的结构，心里都有一张完整的失稳图景。弧长法代码本身并不难，难的是把参数、符号、初始缺陷这些细节一起拿捏住。希望这些经验能帮你少走几次弯路，也希望这篇文章对你有帮助。
本文还有配套的精品资源，点击获取
