Newton–Raphson有限元编程 | 一维非线性弹簧元

知乎、B站：[易木木响叮当]
关注可了解更多的有限元数值仿真技巧。问题或建议，请公众号留言;
如果你觉得木木同学对你有帮助，欢迎赞赏。

今天给大家分享的是：基于Newton–Raphson的一维非线性弹簧有限元分析。

写在前面的话

非线性有限元编程终于有了个开端了，有关这一块要求的功底相对于线弹性编程，着实有点高，不仅要掌握有限元分析的流程，还要熟悉各种数值算法（研一时的数值分析）。不过也不要慌，遇到什么算法就去了解什么算法，千万不要为了掌握有限元中的非线性，而去系统学习数值分析，时间成本太高，除非你有很多的时间(个人理解)。

Kim教授的《Introduction to Nonlinear Finite Element Analysis》是本非常不错的非线性有限元分析教材，难度循序渐进，讲述了有限元分析过程中遇到的各种非线性问题，并提供了相应的Matlab代码，能够让新手小白快速的上手。

书籍地址

Newton–Raphson有限元编程 | 一维非线性弹簧元

问题描述

如下图所示，模型由两个弹簧单元组成，共3个节点，每个节点包含1个自由度，共3个自由度，模型右端受到 的集中载荷，左端固定，弹簧刚度与单元位移量相关，其中， ， ，试求 和 的值。假设单位统一，不计单位。

Newton–Raphson求解

每个弹簧元的平衡方程如下式：

整体刚度组装，去除0左端节点项（固支条件，划行划列）后的整体刚度方程可表示如下：

进一步讲上式化简为1个非线性方程组 ：

公式推导进行到这一步后，Newton–Raphson就可以派上用场了。有限元和数值分析紧密相关的原因就是，有限元思想将系统离散化，通过一系列的理论将物理条件转化为一个个的数学线性/非线性方程，最终依靠数值分析的手段进行问题求解。

有关Newton–Raphson的算法细节，可参考数值分析教材或者木木前几期的推文，本质就是一个切线逼近的过程，本次推文的目的就是带着大家将该算法应用到具体有限元非线性求解案例中。

计算流程

Step1

初值条件：容差 ，最大迭代步=20，初始位移 置零，赋值外载荷 ；

Step2

计算初始残差： ；

收敛参数：

Step3

进入while循环，若 或迭代次数 ，停止循环；

Step4

切线刚度（雅可比矩阵）：

Step5

位移增量 ，参考下图：

Step6

更新位移值： ;

根据新的位移值 ，更新外载荷值 ；

继续更新残差 和收敛参数conv，所用公式可参考Step2；

Step7

更新迭代步：iter = iter + 1，重新进入while循环体中。

代码片

tol = 1.0e-5; iter = 0; 
u = [0; 0];
uold = u;
f = [0; 100];
P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2)
200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)]; % eq. 2.7
R=f-P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2); % eq. 2.12
fprintf('\n iter u1 u2 conv c');
fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv);
while conv > tol && iter < 20
Kt = [600*u(1)+400*u(2)+150 400*(u(1)-u(2))-100
400*(u(1)-u(2))-100 400*u(2)-400*u(1)+100]; % Jacobian matrix
delu = Kt\R; % Fig.2.14
u = uold + delu;
P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2);
200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
R=f -P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
uold = u;
iter = iter + 1;
fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv);
end

结果分析

将计算结果罗列在下表：

由上表结果可知，程序经过5个迭代步最终收敛，

其实该方法已经埋下了一个收敛隐患：解的精度依赖每次迭代求出的雅可比矩阵，而且每次迭代都需要求切线刚度，然后求解线性方程组，使得计算成本较大，为解决该问题，大佬们又在牛顿迭代的基础上进行一些修正，使得计算成本更低，收敛更快、更稳定。该方法可在下次的非线性有限元编程推文中给大家呈现，请大家持续关注~

额外收获

看书的过程中，看到这样一句话：Since stresses or forces are derivatives from displacements in structural problems, it is easy to make the displacement converge rather than the force. 说明了有限元分析过程中，位移控制加载是收敛性比力控制加载好的深层原因，特此记录，分享给大家~

这本书对于非线性非常全面、详细，现在临近毕业，留给我的时间很紧张，只能每次不想写论文的时候，找时间看这本书，度过毕业季后，木木一定将这本书尽可能吃透，然后逐步地分享出来，希望可以帮到有需要的伙伴，最后再一次感谢大家支持木木！