知乎、B站：[易木木響叮當(dāng)]
關(guān)注可了解更多的有限元數(shù)值仿真技巧。問題或建議，請公眾號留言;
如果你覺得木木同學(xué)對你有幫助，歡迎贊賞。

今天給大家分享的主要內(nèi)容：三節(jié)點桿單元Matlab有限元編程，并附帶云圖繪制后處理程序，獲取方法，見文末。

前些日子，木木在學(xué)習(xí)有限元的過程中，在有限元領(lǐng)域發(fā)現(xiàn)一位博主@SAGARBODKHE,他主要也是更新有限元編程的一些視頻，發(fā)現(xiàn)他有一個三節(jié)點桿單元的程序，覺得有趣就分享出來，希望能加深大家對有限元基礎(chǔ)編程的理解，他的主頁：https://www.youtube.com/@SAGARBODKHE，感興趣可以關(guān)注哈~

理論基礎(chǔ)

對比二節(jié)點桿單元

三節(jié)點桿單元二節(jié)點桿單元在單元節(jié)點順序上的區(qū)別是，前者在單元的中心增加一個節(jié)點，節(jié)點順序為：1、3、2。

在位移插值階次的區(qū)別是，前者的位移插值函數(shù)使用的是二次拋物線插值，相比于兩節(jié)點桿單元的線性位移插值函數(shù)。精度更高，可使用較少單元，滿足精度要求。

題外話：

研究生數(shù)值分析教材一般先將誤差分析，然后就是插值擬合，在學(xué)習(xí)插值方法時，也是先講線性插值，然后二次、牛頓等等，在最初學(xué)習(xí)階段，只知插值函數(shù)階次在一定條件下，階次越高，精度越好。

有限元方法與數(shù)值分析思想密切相關(guān)，仔細琢磨一下就可聯(lián)想到位移插值形函數(shù)與數(shù)值分析中拉格朗日插值等思想一致，當(dāng)然還有許多方面也都密切相關(guān)。

位移模式

單元剛度

單元剛度矩陣的推導(dǎo)就不在這里贅述，前人已經(jīng)通過各種方法推導(dǎo)得出剛度矩陣，這里將推導(dǎo)完成的剛度矩陣展示如下：

其中， 為單元長度。

在桿單元的編程中，剛度矩陣不需數(shù)值積分方法求得近似解，直接用現(xiàn)成的即可！

數(shù)值實現(xiàn)

本次的程序沒有設(shè)置函數(shù)文件，直接是一套主程序，運行即可，也可設(shè)置斷點，看每一步如何運行，加上云圖繪制語句共100行左右（包含注釋、空行），行數(shù)較少，可以嘗試去掌握它。

模型描述

位移求解

有限元程序主要有兩部分組成，位移求解、應(yīng)力等其他分量求解及云圖繪制。

雖然云圖繪制并不是必須的，但是效果上很棒，值得學(xué)習(xí)（雖然我不會），但是對已有的繪圖語句會用即可！

Matlab 程序

L=1;
A=1;
E=100e9;
nel=20;

nnp=2*nel+1;
K=zeros(nnp,nnp);
F=zeros(nnp,1);
Lel = L/nel;

%% Stifness Matrix Generation
for i=1:nel
    range=[ 2*i-1 2*i 2*i+1];
    K(range,range)= K(range,range) + (A*E)/(3*Lel)*[7 -8 1;-8 16 -8;1 -8 7];
end

%% Boundary Conditions

fixed_nodes = [1];
free_nodes = setxor(1:nnp,fixed_nodes);%setxor集合運算符（異或）

force_nodes =[nnp];
force_val =[ 1000];
F(force_nodes)= force_val;
%% Partitioning K and F matrix

Kpart = K(free_nodes,free_nodes);
Fpart = F(free_nodes,1);

%% Solve the system of eqn

Up = Kpart\Fpart;
U=zeros(nnp,1);
U(free_nodes)=Up;

在以上程序中需要注意的地方：

1. nnp=2*nel+1該語句中聲明了節(jié)點自由度數(shù)目的大小，大家在套用到2節(jié)點程序時，需要將此改為nnp=nel+1;

2. 在形成總剛時，使用的編程思路和之前我所分享的思路不同，不過我覺得他的這個思路挺方便的，不用另外寫整體剛度的組裝子函數(shù)，大家感興趣可以將該思路與之前我們所采用的整體剛度組裝子程序驗證一下，完全正確；

3. free_nodes = setxor(1:nnp,fixed_nodes)使用到了matlab中setxor集合運算符（異或），兩個集合排除掉中間空白區(qū)域。

其余的編程思路與之前分享推文中的編程思路相同，不加贅述。

位移曲線繪制

%% Visualization
X=zeros(nnp,1);
Y=zeros(nnp,1);

for i=2:nnp
    X(i)=X(i-1) + Lel;
end

Xdef=X+U;
Ydef=Y;
 plot(X,U,'-o')
 title('Displacement variation with length')

云圖繪制

Xdef(end+1)=NaN;
Ydef(end+1)=NaN;
c=U;
c(end+1)=0;

patch(Xdef,Ydef,c,'EdgeColor','interp','Linewidth',5);

colormap jet;
cb = colorbar;
t=get(cb,'Limits');
set(cb,'Ticks',linspace(t(1),t(2),10))
cb.Label.String = 'Resultant Diplacement (m)';
axis equal

xP=zeros(2*nnp,1);
yP=zeros(2*nnp,1);
cP=zeros(2*nnp,1);
h=0.05;
for j=1:2:2*nnp
    xP(j)=Xdef((j+1)/2);
    xP(j+1)=xP(j);
    yP(j)=0+h;
    yP(j+1)=0-h;
    cP(j)=U((j+1)/2);
    cP(j+1)=cP(j);
end

xF=zeros(4,2*nel);
yF=zeros(4,2*nel);
cF=zeros(4,2*nel);
for j=1:2*nel

    xcur=xP((2*j)-1:(2*j)+2);
    ycur=yP((2*j)-1:(2*j)+2);
    ccur=cP((2*j)-1:(2*j)+2);

    xF(:,j)=xcur([1 3 4 2]);
    yF(:,j)=ycur([1 3 4 2]);
    cF(:,j)=ccur([1 3 4 2]);
end

patch(xF,yF,cF)

colormap jet;
cb = colorbar;
t=get(cb,'Limits');
set(cb,'Ticks',linspace(t(1),t(2),10))
cb.Label.String = 'Resultant Diplacement (m)';
axis equal

2 節(jié)點桿單元對比

篇幅原因，這里就不展現(xiàn)2節(jié)點桿單元云圖繪制語句了，可在后臺領(lǐng)取資源文件后，自行嘗試對比。

<<< 左右滑動見更多 >>>

Abaqus對比

為了好玩，我再拿經(jīng)常上手的Abaqus商業(yè)軟件和今天的模型做對比。

<<< 左右滑動見更多 >>>

由上圖可以看到，使用Abaqus對兩節(jié)點、三節(jié)點桿單元進行了有限元模擬，相關(guān)的單元設(shè)置：Mesh-Element type-Truss-Linear/Quadratic。位移值與Matlab計算結(jié)果吻合很好。

Abaqus后處理顯示桿系結(jié)構(gòu)截面形狀

在對桿系結(jié)果進行Abaqus有限元分析時，默認的后處理效果是線條型的。

若想要顯示截面，需要在：View-ODB Display Options-General-Render beam profiles-Scale factor設(shè)置一個看起來合適的系數(shù)即可，如下圖所示：

資源獲取方式

后臺發(fā)送Bar，即可自動獲取源程序。

若覺得對你有用歡迎轉(zhuǎn)發(fā)推文至朋友圈或身邊正在學(xué)習(xí)有限元的小伙伴，感謝大家的支持，點點小贊，在看在看！