材料非线性Matlab有限元编程之初应力法与初应变法（附核心代码）

1年前浏览1900

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

一、切线刚度法

大家可以查阅我上月发布在仿真秀的原创文章《材料非线性Matlab有限元编程之切线刚度法案例应用》这里不再赘述。

二、初应力及初应变法

1、本构方程

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

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

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

2、有限元求解原理

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

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

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

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

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

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

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

3、初应力法

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

（7）

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

（8）

式中，[D]是线弹性材料刚度矩阵，它是非线性材料在时候的切线弹性矩阵。为了使公式7和公式8所表示的应力相同，应当随着的变换，随时调整。比较上述二式，可得

(9)

这里引入假想的线性弹性应力，则

(10)

将公式8代入公式5，得

(11)

令,则

(12)

式中K0就是由线弹性矩阵定义得结构整体刚度矩阵。上式写成矩阵的形式

(13)

(14)

利用上式进行迭代运算，即可求得该非线性问题的解。通常，第一次近似解取为{R0}=0，即使线弹性问题的解。

图4 (a)初应力的含义（b）初应力法迭代过程中单元应力的变化

利用上述方法求解过程中单元中应力应变的变化如图4 所示，最后收敛于真解和。由图中可以看出，如果将当作单元初应力，在图中相当于纵坐标截距。那么整个迭代过程相当于调整所有单元的初应力的过程。即为对应于第n次迭代后初应力场的等效节点力。一旦调整到初应力值时，这个具有初应力场的线弹性解就是非线性弹性问题的解。基于这种物理的解释，此类方法又称为初应力法，或是应力转移法。

因为我们最终目的是要通过有限元编程实现上述方程的求解，所以为了方便编程，将上述初应力迭代方法计算步骤可总结为下述算法步骤

图5 材料非线性问题的初应力迭代算法

4、初应变法

在某些问题中，应力不能由应变显示表达，但应变可以由应力显示表达

(15)

这时，仍采用初应力法会遇到到困难，而采用“初应变”法则比较方便。仿照初应力法，设想用具有初应变的线弹性物理方程来代替公式(15)式表示的非线性物理方程:

(16)

式中为初应变。调整可以使得二式得到相同的应变，比较二式，并令，于是有

(17)

初应变的定义及迭代调整示意图如图6所示。

图6 初应变示意图

将公式16代入公式5，可得

(18)

将上式改写为迭代公式

如果一致位移的第n次近似值为，可以利用公式4求出相应的应变。由及前一次的初应变，利用，可得

再由公式15的本构方程可求得对应的应变值，继而求出新的初应变值

如此迭代多次，直到收敛。初次迭代取，即线弹性解为初始近似值。再迭代过程中单元内应力应变的关系如图6所示。

因为我们最终目的是要通过有限元编程实现上述方程的求解，所以为了方便编程，将上述初应力迭代方法计算步骤可总结为下述算法步骤

图7 材料非线性问题的初应变法迭代算法

4、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, 6) ;       %整体应力矩阵for ie=1:1:element_numberk = 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. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法[bc_number,dummy] = size( gBC1 ) ;for ibc=1:1:bc_numbern = gBC1(ibc, 1 ) ;d = gBC1(ibc, 2 ) ;m = (n-1)*2 + d ;f(m) = gBC1(ibc, 3)* gK(m,m) * 1e20 ;gK(m,m) = gK(m,m) * 1e20 ;end
%% step 5 初应力法迭代E  = gMaterial( gElement(ie, 4), 1 ) ;%生成弹性模量mu = gMaterial( gElement(ie, 4), 2 ) ;%生成泊松比D = [ 1-mu    mu      0mu    1-mu     00      0   (1-2*mu)/2] ;D = D*E/(1-2*mu)/(1+mu) ;             %建立D矩阵gDelta1=ones(node_number * 2,1);      %建立整体 位移向量js=0;Phi1=0;while truefor ie=1:1:element_numberFE = elementforce( ie ,gElementStress) ;AssembleFE( ie, FE ) ;end                                                %组装整体内力向量Phi=( f-gFE );aa=Phi-Phi1;conv=norm(aa)/norm(f)gDelta = gK \ (f - gFE);                           %计算整体 位移向量for ie=1:1:element_numberdelta = NodeDe( ie,gDelta);eps = MatrixB( ie ) * delta;sigma0 = unlinerD(ie,eps)* eps;sigma1 = sigma0 - D * eps;for i = 1:1:3gElementStress(ie,i) = sigma1(i);endend                                                %完成整体到单元的转换，并在单元尺度上计算应变及初应力，并重新生成整体初应力矩阵if abs(max(gDelta1)-max(gDelta))<=1e-1000||js>100               %判断是否达到精度要求breakelsegDelta1=gDelta;Phi1=Phi;endjs=js+1;end

%% step 6. 计算节点应力(采用绕节点加权平均)gNodeStress = zeros( node_number, 6 ) ;       for i=1:node_number                         S = 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
return

上述代码的结构如下：

1. 计算单元刚度矩阵，集成整体刚度矩阵；
2. 计算单元的等效节点力，集成整体节点力向量；
3. 处理约束条件，修改整体刚度矩阵和节点力向量；
4. 初应力法迭代
5. 迭代结束，得到整体节点位移向量；
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

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, 6) ;       %整体应力矩阵for ie=1:1:element_numberk = 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) ) ;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. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法[bc_number,dummy] = size( gBC1 ) ;for ibc=1:1:bc_numbern = gBC1(ibc, 1 ) ;d = gBC1(ibc, 2 ) ;m = (n-1)*2 + d ;f(m) = gBC1(ibc, 3)* gK(m,m) * 1e20 ;gK(m,m) = gK(m,m) * 1e20 ;end%% step 5 初应变法迭代gDelta1=zeros(node_number * 2,1); %取初值delta0=0    E  = gMaterial( gElement(ie, 4), 1 ) ;%生成弹性模量mu = gMaterial( gElement(ie, 4), 2 ) ;%生成泊松比D = [ 1-mu    mu      0mu    1-mu     00      0   (1-2*mu)/2] ;D = D*E/(1-2*mu)/(1+mu) ;             %建立D矩阵%建立整体 位移向量AD=inv(D);%D的逆矩阵js=0;while truefor ie=1:1:element_numberif js==0eps_l=zeros(3,1);enddelta = NodeDe( ie,gDelta1);eps = MatrixB( ie ) * delta; %公式2求epsilon0sigma0= D * (eps-eps_l);epsilon0_sigma=unlinerD_1(ie,sigma0)*sigma0;         %公式3非线性关系epsilon0=epsilon0_sigma-AD*sigma0;eps_l=epsilon0;for i = 1:1:3gElementStrain(ie,i) = epsilon0(i);endFE = elementforce( ie ,gElementStrain(ie,:),D) ;gFE = AssembleFE( ie, FE ) ;endgDelta = gK \ (f + gFE);if abs(max(gDelta1)-max(gDelta))<=1e-1000||js>100               %判断是否达到精度要求breakelsegDelta1=gDelta;endjs=js+1;end   %% step 6. 计算节点应力(采用绕节点加权平均)gNodeStress = zeros( node_number, 6 ) ;       for i=1:node_number                         S = 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) ;endreturn

上述代码的结构如下：

1. 计算单元刚度矩阵，集成整体刚度矩阵；
2. 计算单元的等效节点力，集成整体节点力向量；
3. 处理约束条件，修改整体刚度矩阵和节点力向量；
4. 初应变法迭代
5. 迭代结束，得到整体节点位移向量；
6. 计算单元应力和节点应力；

非线性材料刚度矩阵定义函数如下：

function D_1 = unlinerD_1(ie,sigma)%计算非线性弹性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] ;sigmax = sigma(1);sigmay = sigma(2);D_1 = inv(D)*((E+sigmax+sigmay)/(1-2*mu)/(1+mu) )^(-1);return