材料非线性Matlab有限元编程：初应力法与初应变法

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

一、 切线刚度法

大家可以查阅我上月发布在仿真秀的原创文章《材料非线性Matlab有限元编程：切线刚度法》

二、初应力及初应变法

1、本构方程

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

之所以采用公式1-3三种不同的应力应变关系本构方程，是因为牛顿-拉普森法（切线刚度法）、初应力法、初应变法适用于不同形式的本构方程：切线刚度法，顾名思义，其刚度表达式为应力应变曲线的切线，因此采用微分形式表示其本构关系；初应力法适用于应力由应变确定的本构形式，即应力为应变量，应变为自变量；但某些问题中，应力无法用应变显式表达，相反，应变由应力表达的本构形式，这种情况的非线性本构方程采用初应变法来求解，

图1 材料非线性问题案例模型

2、有限元求解原理

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

图2 三节点三角形单元刚度矩阵推导

但需要注意的是，该平面三角形单元应用的场景是平面应力问题，本案例是平面应变问题，二者的区别如下图所示，除物理方程外，平面应变问题与平面应力问题的变量和方程都完全相同。比较一下这两个物理方程，我们就发现，将平面应力问题里面的弹性模量E换为，把平面应力问题里面的换成，这样的话，我们就从平面应变问题的物理方程就可以转化为平面应变的问题的物理方程，那么反过来也可以由平面应变问题的物理方程换成。因此在对《三角形单元悬臂梁matlab有限元编程》课程代码进行修改的时候，要注意将平面应变问题的材料刚度矩阵，改为平面应变问题的材料刚度矩阵。当然这是针对弹性问题的求解，如果对于材料非线性问题，平面应力应变刚度矩阵是变换的，其本构方程直接采用公式1-3所示的方程来定义。

图3 平面应力问题和平面应变问题的区别

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

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

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

3、初应力法

对于一般非线性材料，物理方程可以表示为

（7）

上式可由具有初应力的线弹性物理方程代替，即