材料非线性Matlab有限元编程之切线刚度法案例应用
导读：《Matlab有限元编程从入门到精通20讲》是笔者原创首发仿真秀的视频课程，本课程为matlab有限元编程专题课，课程主要以案例的形式进行讲解，中间会穿插案例中所涉及到的有限元基本理论，案例不局限于力学问题的有限元求解，还会涉及传热学、电学等问题的有限元求解。因为固体力学领域我最熟悉，所以我们从固体力学开始，所涉及的单元有杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，四面体实体单元等等，力学问题有静力学问题，也有动力学问题，后期还会涉及材料非线性、几何非线性、接触非线性等非线性问题，内容丰富，不断更新完善，预计会有30节课程，每节课都是一个小时左右，干货满满，欢迎订阅！
今天我主要围绕材料非线性问题的有限元Matlab编程求解进行介绍，重点围绕牛顿-拉普森法（切线刚度法）、初应力法、初应变法等三种非线性迭代方法的算法原理展开讲解，最后利用Matlab对材料非线性问题有限元迭代求解算法进行实现，展示了实现求解的核心代码。由于篇幅原因，本文分两期分享，今天重点介绍切线刚度法（牛顿-拉普森法）案例应用。
一、切线刚度法（牛顿-拉普森法）案例背景
本次课程同样以具体案例为对象进行材料非线性问题的有限元变成求解，求解模型如图1，模型边界为20m×10m，材料本构方程如公式1所示，其中，弹性模量E=20MPa，泊松比0.35，模型上表面中间位置作用20kPa超载，超载作用范围为4m。按照平面应变问题考虑，使用常应变三角形单元分析模型上表面中间点竖向沉降，对应的有限元模型和计算结果如图2、3所示。
（1）
图1 两物体接触的平面应力问题
图2 有限元模型
图3 切线刚度法计算结果
二、有限元求解原理
由图2所示，有限元离散方式采用的是三节点三角形单元进行离散，因此我们要有三角形平面单元弹性问题的求解基础知识，大家可以观看仿真秀视频课程《Matlab有限元编程从入门到精通20讲》课程中的“三角形单元悬臂梁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 SolveModelglobal 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_numberenf = 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=0gDelta=zeros(node_number * 2,1); %取初值delta0=0js=0;while truegK=zeros( node_number * 2, node_number * 2 );gFE=zeros(node_number * 2,1);for ie=1:1:element_numberdelta = NodeDe( ie,gDelta1);eps = MatrixB( ie ) * delta; %公式2求epsilon0sigma0 = unlinerD(ie,eps)* (eps-gElementStrain(ie,:)')+gElementStress(ie,:)';%g公式3求sigmma0kt = 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:3gElementStrain(ie,ii) = eps(ii);gElementStress(ie,ii) = sigma0(ii);endFE = elementforce( ie ,gElementStress) ;AssembleFE( ie, FE ) ;end;Phi=( f-gFE );%% step 5. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法[bc_number,dummy] = size( gBC1 ) ;for ibc=1:1:bc_numbern = 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 ;enddd = gK \ Phi;gDelta=dd+gDelta;conv=norm(Phi)/norm(f)if js>100 || conv<1e-4 %判断是否达到精度要求breakelsegDelta1=gDelta;endend%% step 6. 计算节点应力(采用绕节点加权平均)gNodeStress = zeros( node_number, 6 ) ;for i=1:node_numberS = zeros( 1, 3 ) ;A = 0 ;for ie=1:1:element_numberfor k=1:1:3if i == gElement( ie, k )area= ElementArea( ie ) ;S = S + gElementStress(ie,1:3 ) * area ;A = A + area ;break ;endendendgNodeStress(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
6、计算单元应力和节点应力；
非线性材料刚度矩阵定义函数如下：
function D = unlinerD (ie,eps)%计算非线性弹性D矩阵% 输入参数:% ie ---- 单元号% 返回值:% D ---- D矩阵global gElement gMaterialE = gMaterial( gElement(ie, 4), 1 ) ;mu = gMaterial( gElement(ie, 4), 2 ) ;D = [ 1-mu mu 0mu 1-mu 00 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有限元编程从入门到精通20讲》后面进行加餐，欢迎识别下方二维码试看。
我的Matlab有限元编程精品课
Matlab有限元编程从入门到精通20讲：快速获得各典型有限元案例的Matlab代码
可开电子发票，赠送答疑专栏
提供vip群交流，课程可反复回看
本课程为matlab有限元编程专题课，课程主要以案例的形式进行讲解，中间会穿插案例中所涉及到的有限元基本理论，案例不局限于力学问题的有限元求解，还会涉及传热学、电学等问题的有限元求解。
因为固体力学领域我最熟悉，所以我们从固体力学开始，所涉及的单元有杆单元，梁单元，平面三角形单元，薄板单元，厚板单元，四面体实体单元等等，力学问题有静力学问题，也有动力学问题，后期还会涉及材料非线性、几何非线性、接触非线性等非线性问题，内容丰富，不断更新完善。
此外，笔者为所有订阅用户提供知识圈答疑服务和VIP用户交流群。并附赠课程相关资料等（平台支持自行开具电子发票）。
参考文献：
