材料非线性Matlab有限元编程：切线刚度法（牛顿-拉普森法）

导读：

 本文主要围绕材料非线性问题的有限元Matlab编程求解进行介绍，重点围绕牛顿-拉普森法（切线刚度法）、初应力法、初应变法等三种非线性迭代方法的算法原理展开讲解，最后利用Matlab对材料非线性问题有限元迭代求解算法进行实现，展示了实现求解的核心代码。这些内容都将收录在我的原创精品课《matlab有限元编程从入门到精通》

【课程完整案例源码链接】：接触有限元编程 | Matlab源码 | 理论文本

今天我主要围绕材料非线性问题的有限元Matlab编程求解进行介绍，重点围绕牛顿-拉普森法（切线刚度法）、初应力法、初应变法等三种非线性迭代方法的算法原理展开讲解，最后利用Matlab对材料非线性问题有限元迭代求解算法进行实现，展示了实现求解的核心代码。由于篇幅原因，本文分两期分享，今天重点介绍切线刚度法（牛顿-拉普森法）案例应用。

一、切线刚度法（牛顿-拉普森法）案例背景

本次课程同样以具体案例为对象进行材料非线性问题的有限元变成求解，求解模型如图1，模型边界为20m×10m，材料本构方程如公式1所示，其中，弹性模量E=20MPa，泊松比0.35，模型上表面中间位置作用20kPa超载，超载作用范围为4m。按照平面应变问题考虑，使用常应变三角形单元分析模型上表面中间点竖向沉降，对应的有限元模型和计算结果如图2、3所示。

（1）

图1 两物体接触的平面应力问题

图2 有限元模型

图3 切线刚度法计算结果

二、有限元求解原理

由图2所示，有限元离散方式采用的是三节点三角形单元进行离散，因此我们要有三角形平面单元弹性问题的求解基础知识，大家可以观看仿真秀视频课程《matlab有限元编程从入门到精通》中的“三角形单元悬臂梁matlab有限元编程”小节，详细讲解了基于三角形三节点单元的有限元离散过程以及弹性刚度矩阵的推导

图4 三角形单元悬臂梁matlab有限元编程截图

在掌握基于三角形单元弹性问题的求解基础知识后，针对本案例的纯材料非线性问题，其几何方程、平衡方程的建立均为线性关系，只有物理方程存在非线性关系，具体分析如下：属于小变形问题，因此公式2表示的几何关系是线性的，公式3以应力形式表示的平衡条件也是线性的。引入物理方程，其一般形式为：

（2）

（3）

在材料非线性问题中，应力与应变关系是非线性的，对于本案例，应力应变的关系如公式1所示。所以，以节点位移列阵表示的平衡方程不再是线性的，可以写成

（4）

上式与几何非线性的的表达式类似，因此材料非线性和几何非线性都可以用相同的迭代方法来求解。本系列课程主要介绍牛顿-拉普森法（切线刚度法）、初应力法、初应变法等三种迭代方法。这一小节围绕最常用的非线性迭代算法——切线刚度发（牛顿拉普森法）进行介绍。

三、牛顿-拉普森法（切线刚度法）

如果材料的应力应变关系能表示成公式5所示的增量形式，就可以利用切线刚度法。

（5）

式中的D_t表示材料的切线弹性矩阵，改写公式3可得

（6）

由于式中B和R均为定常量，对上式进行微分计算可得，

（7）

将公式5和公式2带入公式7可得

（8）

（9）

上式K_T为系统的切线刚度矩阵。

利用牛顿-拉普森方法（切线刚度法），则迭代公式为

（10）

式中为第n次近似解时，等效节点里和等效节点荷载的差，称为失衡力。当失衡力在精度范围内时，则认为迭代收敛，获得了方程式6的解。

上述牛顿-拉普森迭代方法计算步骤可总结为下图

图5 材料非线性问题的牛顿-拉普森迭代算法

四、Matlab程序设计

这里我们展示了求解该有限元模型的核心代码，主要涉及牛顿-拉普森迭代算法以及非线性材料刚度矩阵的定义。

function SolveModel
global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gDF gElementStrain gFE
%% step 1. 定义整体刚度矩阵和节点力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 2, node_number * 2 ) ;
gFE = sparse( node_number * 2, 1 ) ;               %整体内力向量
f = sparse( node_number * 2, 1 ) ;    
%% step 2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
[element_number, dunmmy] = size( gElement ) ;
gElementStrain = zeros( element_number, 3) ;       %整体应变矩阵
gElementStress = zeros( element_number, 3) ;       %整体应力矩阵
% for ie=1:1:element_number
% k = StiffnessMatrix( ie ) ;
% AssembleStiffnessMatrix( ie, k ) ;
% end
%% step 3. 计算地面超载产生的等效节点力
[df_number,dummy] = size( gDF ) ;
for idf = 1:1:df_number
enf = EquivalentNodeForce( gDF(idf,1), gDF(idf,2), gDF(idf,3), gDF(idf,4) )*10;
i = gElement( gDF(idf,1), 1 ) ;
j = gElement( gDF(idf,1), 2 ) ;
m = gElement( gDF(idf,1), 3 ) ;
f( (i-1)*2+1 : (i-1)*2+2 ) = f( (i-1)*2+1 : (i-1)*2+2 ) + enf( 1:2 ) ;
f( (j-1)*2+1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + enf( 3:4 ) ;
f( (m-1)*2+1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + enf( 5:6 ) ;
end         
%% step 4 切线刚度法迭代
gDelta1=zeros(node_number * 2,1); %取初值delta0=0    
gDelta=zeros(node_number * 2,1); %取初值delta0=0    
js=0;
while true
gK=zeros( node_number * 2, node_number * 2 );
gFE=zeros(node_number * 2,1);
for ie=1:1:element_number        
delta = NodeDe( ie,gDelta1);
eps = MatrixB( ie ) * delta; %公式2求epsilon0                      
sigma0 = unlinerD(ie,eps)* (eps-gElementStrain(ie,:)')+gElementStress(ie,:)';%g公式3求sigmma0
kt = zeros( 6, 6 ) ;
h  = gMaterial( gElement(ie, 4), 3 ) ;
B = MatrixB( ie );
area = ElementArea( ie );
kt = transpose(B)*unlinerD(ie,eps)*B*h*abs(area) ;    
AssembleStiffnessMatrix( ie, kt ) ;
for ii = 1:1:3
gElementStrain(ie,ii) = eps(ii);
gElementStress(ie,ii) = sigma0(ii);
end
FE = elementforce( ie ,gElementStress) ;
AssembleFE( ie, FE ) ; 
end;
Phi=( f-gFE );
%% step 5. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*2 + d ;
Phi(m) = gBC1(ibc, 3)* gK(m,m) * 1e20 ;
gK(m,m) = gK(m,m) * 1e20 ;
end
dd = gK \ Phi;  
gDelta=dd+gDelta;
conv=norm(Phi)/norm(f)
if js>100 || conv<1e-4              %判断是否达到精度要求
break
else
gDelta1=gDelta;     
end
end
%% step 6. 计算节点应力(采用绕节点加权平均)
gNodeStress = zeros( node_number, 6 ) ;       
for i=1:node_number                         
S = zeros( 1, 3 ) ;                         
A = 0 ;
for ie=1:1:element_number
for k=1:1:3
if i == gElement( ie, k ) 
area= ElementArea( ie ) ;
S = S + gElementStress(ie,1:3 ) * area ;
A = A + area ;
break ;
end
end
end
gNodeStress(i,1:3) = S / A ;
gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
end

上述代码的结构如下：1、计算单元刚度矩阵，集成整体刚度矩阵；2、计算单元的等效节点力，集成整体节点力向量；3、处理约束条件，修改整体刚度矩阵和节点力向量；4、切线刚度法迭代5、迭代结束，得到整体节点位移向量；

其中上述计算单元力的应力是需要的非线性材料刚度矩阵定义函数如下：

function D = unlinerD (ie,eps)
%计算非线性弹性D矩阵
% 输入参数:
% ie ----  单元号
% 返回值:
% D  ----  D矩阵
global  gElement gMaterial
E  = gMaterial( gElement(ie, 4), 1 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
D = [ 1-mu    mu      0
mu    1-mu     0
0      0   (1-2*mu)/2] ;
epsx = eps(1);
epsy = eps(2);
D = D*E*(1-100*epsx^2-100*epsy^2)/(1-2*mu)/(1+mu) ;
return

五、Matlab程序组成

以上就是笔者围绕牛顿-拉普森法（切线刚度法）进行的讲解，接下来还会介绍初应力法、初应变法等二种非线性迭代方法的算法原理展开讲解，欢迎大家持续关注。由于篇幅原因至列举部分源码，真正实现上述问题的求解还需要包括函数具体如下，整个项目的求解源码发布在《Matlab有限元编程从入门到精通30讲》课程资料附件中。需要的小伙伴也可以私信我~

我的Matlab有限元编程精品课

本课程为matlab有限元编程专题课，课程主要以案例的形式进行讲解，中间会穿插案例中所涉及到的有限元基本理论，案例不局限于力学问题的有限元求解，还会涉及传热学、电学等问题的有限元求解。

因为固体力学领域我最熟悉，所以我们从固体力学开始，所涉及的单元有杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，四面体实体单元等等，力学问题有静力学问题，也有动力学问题，后期还会涉及材料非线性、几何非线性、接触非线性等非线性问题，内容丰富，不断更新完善。