非线性有限元编程 | 增量迭代法
在介绍了前几期推文后，相信大家对非线性算法有着一定的了解，各有优缺点及适用情景，但毫无疑问的是在保证初值选定合适的情况下，牛顿法收敛速度最快。
在固体力学问题中，初始估计通常设定为结构的未变形形状，也就是说，假定所有位移最初都是零。由于初始估计通常从零位移开始，当施加载荷较小时，Newton - Raphson方法快速收敛到解，收敛困难发生在大的应用载荷，造成大的位移幅度。
于是一个idea就此产生，那我们在系统中施加载荷时，分段施加，每一段施加相对较小的载荷，在每个载荷段内，Newton - Raphson方法都可以获得收敛解。
该方法也就是我们在商软中常常见到的增量迭代法，本期推文将分别从力控制加载和位移控制加载两个角度来解析增量迭代法。
【声明】：本次案例分享来自Kim教授的《Introduction to Nonlinear Finite Element Analysis》。
力控制加载
考虑一单变量问题，如下图所示，在第一个增量中，假设所施加的荷载为 ，从初始估计的 开始求解非线性方程，选择这个增量的大小是为了使数值方法能够快速收敛到解 。在第二个增量中，所应用的负载增加到 ，初始位移更新为 ，同样，选择增量的大小 是为了使数值方法能够快速收敛到解 。重复上述过程，直到所应用的载荷增量达到总载荷值。
注：在分段施加载荷时，没必要平均施加。
位移控制加载
在外界荷载作用下，材料内部出现软化行为（弹塑性材料），如下图所示，使用力控制加载法很难捕捉到软化曲线段，若施加的载荷大于 ，那么求解迭代过程将不会收敛。
这时有了新的idea，纵坐标分段不行，那就横坐标分段施加，具体怎么实现的呢？
力控制加载时，求解出的物理量是节点位移；当位移控制加载时（分段施加位移），求解出的物理量是节点力。在数学上，这两个过程是等价的。
但实际上，位移控制的过程比力控制的过程更稳定。原因就如上图所示，当材料内出现软化行为时，力法显得有点力不从心，位移法却能迎刃而解，扩大了收敛解的范围，当然介绍的只是其中一种情况，位移控制加载的优势相信大家在使用商软的过程中，已经体会到了。
接下来，将位移控制加载法应用于非线性弹簧案例中，与之前的力控制加载结果做对比。
FEM非线性案例
如下图所示，模型由两个弹簧单元组成，共3个节点，每个节点包含1个自由度，共3个自由度，模型右端受到 的集中载荷，左端固定，弹簧刚度与单元位移量相关，其中， ， ，试求 和 的值。假设单位统一，不计单位。
每个弹簧元的平衡方程如下式：
整体刚度组装，去除0左端节点项（固支条件，划行划列）后的整体刚度方程可表示如下：
进一步讲上式化简为1个非线性方程组 ：
求解代码
clear
tol = 1.0e-5; conv = 0; u1 = 0; u1old = u1;
fprintf('\n step u1 u2 F');
uu1 = zeros(9,1);
uu2 = zeros(9,1);
FF = zeros(9,1);
% Displacement increment loop
for i=1:9
u2 = 0.1*i;
uu2(i) = u2;
P = 300*u1^2+400*u1*u2-200*u2^2+150*u1-100*u2;
R = -P;
conv = R^2;
% Convergence loop
iter = 0;
while conv > tol && iter < 20
Kt = 600*u1+400*u2+150;
delu1 = R/Kt;
u1 = u1old + delu1;
uu1(i) = u1;
P = 300*u1^2+400*u1*u2-200*u2^2+150*u1-100*u2;
R = -P;
conv= R^2;
u1old = u1;
iter = iter + 1;
end
F = 200*u1^2-400*u1*u2+200*u2^2-100*u1+100*u2;
FF(i) = F;
fprintf('\n %3d %7.5f %7.5f %7.3f',i,u1,u2,F);
end
plot(FF,uu1,'-O',FF,uu2,'-s')
在以上程序中，设置了位移以0.1为增量的形式，从0.1增加至0.9，最终求解得到的节点力为100，与原载荷吻合，而且在求解过程中，每一步均发生了一次迭代，极为迅速。
思考
在介绍了今天的推文后，有没有体会到有限元求解过程中，最开始我们是学会了如何怎样求解，有了一定了解后，我们又在想面对非线性问题如何求解的更快，诞生了各种非线性算法，在算法的基础上产生新想法：增量迭代，后续在面对大规模计算时，为加快求解速度，可能还会接触到并行计算。学习过程亦是如此，由简入难，随着你学习程度的加深，再回过头看当初觉得很难的问题，可能会觉得“好简单”。
