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