導讀：本文主要圍繞材料非線性問題的有限元Matlab編程求解進行介紹，重點圍繞牛頓-拉普森法（切線剛度法）、初應力法、初應變法等三種非線性迭代方法的算法原理展開講解，最后利用Matlab對材料非線性問題有限元迭代求解算法進行實現(xiàn)，展示了實現(xiàn)求解的核心代碼。這些內(nèi)容都將收錄在我的原創(chuàng)精品課《matlab有限元編程從入門到精通》。
【課程完整案例源碼鏈接】：材料非線性有限元編程|切線剛度法 | 初應力法 | 初應變法【Matlab源碼+理論文本】
一、 切線剛度法
大家可以查閱我上月發(fā)布在仿真秀的原創(chuàng)文章《材料非線性Matlab有限元編程：切線剛度法》
二、初應力及初應變法
1、本構方程
本文同樣以具體案例為對象進行材料非線性問題的有限元變成求解，求解模型如圖1，模型邊界為20m×10m，公式1-3材料本構方程如公式1所示，其中，彈性模量E=20MPa，泊松比0.35，模型上表面中間位置作用20kPa超載，超載作用范圍為4m。按照平面應變問題考慮，使用常應變?nèi)切螁卧治瞿Ｐ蜕媳砻嬷虚g點豎向沉降，對應的有限元模型和計算結果如圖2、3所示。
之所以采用公式1-3三種不同的應力應變關系本構方程，是因為牛頓-拉普森法（切線剛度法）、初應力法、初應變法適用于不同形式的本構方程：切線剛度法，顧名思義，其剛度表達式為應力應變曲線的切線，因此采用微分形式表示其本構關系；初應力法適用于應力由應變確定的本構形式，即應力為應變量，應變?yōu)樽宰兞浚坏承﹩栴}中，應力無法用應變顯式表達，相反，應變由應力表達的本構形式，這種情況的非線性本構方程采用初應變法來求解，
圖1 材料非線性問題案例模型
2、有限元求解原理
由圖2所示，有限元離散方式采用的是三節(jié)點三角形單元進行離散，因此我們要有三角形平面單元彈性問題的求解基礎知識，大家可以觀看b站的《Matlab有限元編程從入門到精通》課程中的“三角形單元懸臂梁matlab有限元編程”小節(jié)，詳細講解了基于三角形三節(jié)點單元的有限元離散過程以及彈性剛度矩陣的推導。
圖2 三節(jié)點三角形單元剛度矩陣推導
但需要注意的是，該平面三角形單元應用的場景是平面應力問題，本案例是平面應變問題，二者的區(qū)別如下圖所示，除物理方程外，平面應變問題與平面應力問題的變量和方程都完全相同。比較一下這兩個物理方程，我們就發(fā)現(xiàn)，將平面應力問題里面的彈性模量E換為，把平面應力問題里面的換成，這樣的話，我們就從平面應變問題的物理方程就可以轉化為平面應變的問題的物理方程，那么反過來也可以由平面應變問題的物理方程換成。因此在對《三角形單元懸臂梁matlab有限元編程》課程代碼進行修改的時候，要注意將平面應變問題的材料剛度矩陣，改為平面應變問題的材料剛度矩陣。當然這是針對彈性問題的求解，如果對于材料非線性問題，平面應力應變剛度矩陣是變換的，其本構方程直接采用公式1-3所示的方程來定義。
圖3 平面應力問題和平面應變問題的區(qū)別
在掌握基于三角形單元彈性問題的求解基礎知識后，針對本案例的純材料非線性問題，其幾何方程、平衡方程的建立均為線性關系，只有物理方程存在非線性關系，具體分析如下：屬于小變形問題，因此公式2表示的幾何關系是線性的，公式3以應力形式表示的平衡條件也是線性的。引入物理方程，其一般形式為
在材料非線性問題中，應力與應變關系是非線性的，對于本案例，應力應變的關系如公式1所示。所以，以節(jié)點位移列陣表示的平衡方程不再是線性的，可以寫成
上式與幾何非線性的的表達式類似，因此材料非線性和幾何非線性都可以用相同的迭代方法來求解。本系列課程主要介紹牛頓-拉普森法（切線剛度法）、初應力法、初應變法等三種迭代方法。這一小節(jié)圍繞初應力和初應變法進行介紹。
3、初應力法
對于一般非線性材料，物理方程可以表示為
（7）
上式可由具有初應力的線彈性物理方程代替，即
添加圖片注釋，不超過 140 字（可選）
（8）
式中，[D]是線彈性材料剛度矩陣，它是非線性材料在 時候的切線彈性矩陣。為了使公式7和公式8所表示的應力相同，應當隨著 的變換，隨時調(diào)整 。比較上述二式，可得
添加圖片注釋，不超過 140 字（可選）
(9)
這里引入假想的線性彈性應力
添加圖片注釋，不超過 140 字（可選）
，則
添加圖片注釋，不超過 140 字（可選）
(10)
將公式8代入公式5，得
添加圖片注釋，不超過 140 字（可選）
(11)
令
添加圖片注釋，不超過 140 字（可選）
,則
添加圖片注釋，不超過 140 字（可選）
(12)
式中K0就是由線彈性矩陣定義得結構整體剛度矩陣。上式寫成矩陣的形式
添加圖片注釋，不超過 140 字（可選）
(13)
添加圖片注釋，不超過 140 字（可選）
(14)
利用上式進行迭代運算，即可求得該非線性問題的解。通常，第一次近似解取為{R0}=0，即使線彈性問題的解。
圖4 (a)初應力的含義 （b）初應力法迭代過程中單元應力的變化
利用上述方法求解過程中單元中應力應變的變化如圖4 所示，最后收斂于真解 和 。由圖中可以看出，如果將
添加圖片注釋，不超過 140 字（可選）
當作單元初應力，在圖中相當于縱坐標截距。那么整個迭代過程相當于調(diào)整所有單元的初應力的過程。 即為對應于第n次迭代后初應力場 的等效節(jié)點力。一旦調(diào)整到初應力值
添加圖片注釋，不超過 140 字（可選）
時，這個具有初應力場的線彈性解就是非線性彈性問題的解。基于這種物理的解釋，此類方法又稱為初應力法，或是應力轉移法。
因為我們最終目的是要通過有限元編程實現(xiàn)上述方程的求解，所以為了方便編程，將上述初應力迭代方法計算步驟可總結為下述算法步驟
圖5 材料非線性問題的初應力迭代算法
4、初應變法
在某些問題中，應力不能由應變顯示表達，但應變可以由應力顯示表達
添加圖片注釋，不超過 140 字（可選）
(15)
這時，仍采用初應力法會遇到到困難，而采用“初應變”法則比較方便。仿照初應力法，設想用具有初應變的線彈性物理方程來代替公式(15)式表示的非線性物理方程:
(16)
式中 為初應變。調(diào)整可以使得二式得到相同的應變，比較二式，并令
，于是有
(17)
初應變的定義及迭代調(diào)整示意圖如圖6所示。
圖6 初應變示意圖
將公式16代入公式5，可得
(18)
將上式改寫為迭代公式
（19）
如果一致位移的第n次近似值為 ，可以利用公式4求出相應的應變 。由 及前一次的初應變 ，利用
（20）
，可得
（21）
再由公式15的本構方程可求得對應 的應變值 ，繼而求出新的初應變值
(22)
如此迭代多次，直到收斂。初次迭代取
(23)
，即線彈性解為初始近似值。再迭代過程中單元內(nèi)應力應變的關系如圖6所示。
因為我們最終目的是要通過有限元編程實現(xiàn)上述方程的求解，所以為了方便編程，將上述初應力迭代方法計算步驟可總結為下述算法步驟
圖7 材料非線性問題的初應變法迭代算法
4、Matlab程序設計-初應力法
這里我們展示了求解該有限元模型的核心代碼，主要涉及初應力法和初應變法的迭代算法以及非線性材料剛度矩陣的定義。
下述代碼為初應力法的迭代算法的實現(xiàn)。
function SolveModel
global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gDF gElementStrain gFE
%% step 1. 定義整體剛度矩陣和節(jié)點力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 2, node_number * 2 ) ;
gFE = sparse( node_number * 2, 1 ) ; %整體內(nèi)力向量
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_number
k = StiffnessMatrix( ie ) ;
AssembleStiffnessMatrix( ie, k ) ;
end
%% step 3. 計算地面超載產(chǎn)生的等效節(jié)點力
[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. 處理約束條件，修改剛度矩陣和節(jié)點力向量。采用乘大數(shù)法
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
n = 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 0
mu 1-mu 0
0 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 true
for ie=1:1:element_number
FE = elementforce( ie ,gElementStress) ;
AssembleFE( ie, FE ) ;
end %組裝整體內(nèi)力向量
Phi=( f-gFE );
aa=Phi-Phi1;
conv=norm(aa)/norm(f)
gDelta = gK \ (f - gFE); %計算整體位移向量
for ie=1:1:element_number
delta = NodeDe( ie,gDelta);
eps = MatrixB( ie ) * delta;
sigma0 = unlinerD(ie,eps)* eps;
sigma1 = sigma0 - D * eps;
for i = 1:1:3
gElementStress(ie,i) = sigma1(i);
end
end %完成整體到單元的轉換，并在單元尺度上計算應變及初應力，并重新生成整體初應力矩陣
if abs(max(gDelta1)-max(gDelta))<=1e-1000||js>100 %判斷是否達到精度要求
break
else
gDelta1=gDelta;
Phi1=Phi;
end
js=js+1;
end
%% step 6. 計算節(jié)點應力(采用繞節(jié)點加權平均)
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
return上述代碼的結構如下：
非線性材料剛度矩陣定義函數(shù)如下：
function D = unlinerD (ie,eps)
%計算非線性彈性D矩陣
% 輸入?yún)?shù):
% 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) ;
return5、Matlab程序設計-初應力法
這里我們展示了求解該有限元模型的核心代碼，主要涉及初應力法和初應變法的迭代算法以及非線性材料剛度矩陣的定義。下述代碼為初應力法的迭代算法的實現(xiàn)。
function SolveModel
global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gDF gElementStrain gFE
%% step 1. 定義整體剛度矩陣和節(jié)點力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 2, node_number * 2 ) ;
gFE = sparse( node_number * 2, 1 ) ; %整體內(nèi)力向量
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_number
k = StiffnessMatrix( ie ) ;
AssembleStiffnessMatrix( ie, k ) ;
end
%% step 3. 計算地面超載產(chǎn)生的等效節(jié)點力
[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) ) ;
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. 處理約束條件，修改剛度矩陣和節(jié)點力向量。采用乘大數(shù)法
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
n = 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 0
mu 1-mu 0
0 0 (1-2*mu)/2] ;
D = D*E/(1-2*mu)/(1+mu) ; %建立D矩陣
%建立整體位移向量
AD=inv(D);%D的逆矩陣
js=0;
while true
for ie=1:1:element_number
if js==0
eps_l=zeros(3,1);
end
delta = NodeDe( ie,gDelta1);
eps = MatrixB( ie ) * delta; %公式2求epsilon0
sigma0= D * (eps-eps_l);
epsilon0_sigma=unlinerD_1(ie,sigma0)*sigma0; %公式3非線性關系
epsilon0=epsilon0_sigma-AD*sigma0;
eps_l=epsilon0;
for i = 1:1:3
gElementStrain(ie,i) = epsilon0(i);
end
FE = elementforce( ie ,gElementStrain(ie,:),D) ;
gFE = AssembleFE( ie, FE ) ;
end
gDelta = gK \ (f + gFE);
if abs(max(gDelta1)-max(gDelta))<=1e-1000||js>100 %判斷是否達到精度要求
break
else
gDelta1=gDelta;
end
js=js+1;
end
%% step 6. 計算節(jié)點應力(采用繞節(jié)點加權平均)
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
return上述代碼的結構如下：
非線性材料剛度矩陣定義函數(shù)如下：
function D_1 = unlinerD_1(ie,sigma)
%計算非線性彈性D矩陣
% 輸入?yún)?shù):
% 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] ;
sigmax = sigma(1);
sigmay = sigma(2);
D_1 = inv(D)*((E+sigmax+sigmay)/(1-2*mu)/(1+mu) )^(-1);
return三、我的Matlab有限元編程視頻教程
以上就是筆者圍繞材料非線性Matlab有限元編程之初應力法與初應變法進行的講解，
由于篇幅原因至列舉部分源碼，真正實現(xiàn)上述問題的求解還需要包括函數(shù)具體如下，整個項目的求解源碼發(fā)布在《Matlab有限元編程從入門到精通30講》課程資料附件中。需要的小伙伴也可以私信我~
我的Matlab有限元編程精品課
本課程為matlab有限元編程專題課，課程主要以案例的形式進行講解，中間會穿插案例中所涉及到的有限元基本理論，案例不局限于力學問題的有限元求解，還會涉及傳熱學、電學等問題的有限元求解。
因為固體力學領域我最熟悉，所以我們從固體力學開始，所涉及的單元有桿單元，梁單元，平面三角形單元，薄板單元，厚板單元，四面體實體單元等等，力學問題有靜力學問題，也有動力學問題，后期還會涉及材料非線性、幾何非線性、接觸非線性等非線性問題，內(nèi)容豐富，不斷更新完善。