非线性有限元编程 | Modified Newton–Raphson

本次分享的是基于Modified Newton–Raphson的非线性有限元编程小案例，将会介绍牛顿法的优缺点，以及如何一步步修正，最后结合弹簧单元的例子进行Modified Newton–Raphson数值实现，望大家喜欢。

牛顿法缺点

考虑使用N-R法求解以下非线性方程：

该方程的精确解为 ，设置收敛容差 ，初始值 。

每次迭代的切线：

牛顿法求解代码

xdata=zeros(40,1);
ydata=zeros(40,1);
tol = 1.0e-5;
iter = 0;
u = 0.5;
uold = u;
c=0;
P = u+atan(5*u);
R = -P;
conv= R^2;
xdata(1)=u;
ydata(1)=P;
while conv > tol && iter < 20
Kt = 1+5*(cos(atan(5*u)))^2;
delu = R/Kt;
u = uold + delu;
P = u+atan(5*u);
R = -P;
conv= R^2;
uold = u;
iter = iter + 1;
xdata(2*iter)=u; ydata(2*iter)=0;
xdata(2*iter+1)=u; ydata(2*iter+1)=P;
end
%
plot(xdata,ydata);
hold on;
x=[-1:0.1:1];
y=x+atan(5*x);
plot(x,y)

以上程序在迭代次数超过预置的20次后停止迭代，每次迭代的结果皆偏离精确解，导致数值发散。

不收敛原因：初始取值不合理

解决方法：初值小于0.5，程序将会很快收敛；大于0.5，程序发散。

由此可以得出Newton–Raphson method的优缺点：

1. 收敛快
2. 每次都要更新切线刚度，计算成本高
3. 初始值在接近精确值时。才能保证收敛

修正牛顿法

为克服牛顿法在求解时的不足，可以将牛顿法的求解思路稍加改进，本节将简要回顾一下Modified Newton–Raphson method的思想。相比于Newton–Raphson method，Modified的区别：

1. 每次迭代只需 the initial  tangent stiffness matrix
2. 迭代次数多，每次迭代时间长，但计算成本相对较低

为了提高收敛性，可先迭代一些tangent stiffness matrix，然后进行Modified Newton–Raphson method。但该方法在构建合理的tangent stiffness matrix，很难指定确定迭代次数。

非线性应用

Example

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

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

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

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

Solution

tol = 1.0e-5;
iter = 0;
u = [0.2; 0.4];
uold = u;
c=0;
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)];
R=f-P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
fprintf('\n iter u1 u2 conv c');
fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
Kt = [600*u(1)+400*u(2)+150 -400*u(2)+400*u(1)-100
400*u(1)-400*u(2)-100 400*u(2)-400*u(1)+100];
while conv > tol && iter < 100
delu = Kt\R;
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);
c = abs(0.9-u(2))/abs(0.9-uold(2))^2;
uold = u;
iter = iter + 1;
fprintf('\n %3d %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end

由以上求解过程，可以看出，Modified Newton–Raphson method的迭代次数相对于Newton–Raphson method的迭代次数较多，问题规模较为简单，体现不出求解时间的差别，读者可自行验证。

【注】:在本程序中，设置的初值是u = [0.2; 0.4]，读者可设置不同的初值，探究初值对收敛次数的影响。在设置初值时，应保持 (常识判断)，若反向设置，将导致程序错误。