如何用matlab實(shí)現(xiàn)針對(duì)四面體單元?jiǎng)澐值娜S結(jié)構(gòu)進(jìn)行有限元編程

導(dǎo)讀：自9月以來(lái)，我的《教你Matlab有限元編程對(duì)懸臂梁進(jìn)行受力分析》原創(chuàng)文章發(fā)布以來(lái)，我的精品課《Matlab有限元編程從入門(mén)到精通10講》受到了越來(lái)越多的工程師朋友的關(guān)注，已超過(guò)50多人已經(jīng)訂閱學(xué)習(xí)，歡迎大家加入訂閱用戶交流群一起討論Matlab有限元編程那些事（哈哈哈，為大家答疑和分享資料）。

今天，我接著分享Matlab有限元編程專業(yè)技能。之前的課程我們學(xué)習(xí)了一維梁?jiǎn)卧?，二維平面單元，三維板殼單元的matlab有限元編程，本次案例主要講解如何用matlab實(shí)現(xiàn)針對(duì)四面體單元?jiǎng)澐值娜S結(jié)構(gòu)進(jìn)行有限元編程，具體案例是一個(gè)懸臂梁受集中荷載的問(wèn)題。圖1為本案例Matlab編程計(jì)算得到的結(jié)果。主要內(nèi)容涉及四面體單元的有限元基本理論的推導(dǎo)，主要是單元?jiǎng)偠染仃嚨耐茖?dǎo)，此外還包括等參單元和Hammer數(shù)值積分以及三維問(wèn)題的后處理計(jì)算。

圖1 懸臂梁受集中荷載的應(yīng)力云圖

一個(gè)完整的有限元程序基本組成部分包括前處理模塊、分析主程序模塊和后處理模塊。在前處理模塊中，實(shí)現(xiàn)節(jié)點(diǎn)坐標(biāo)輸入、單元節(jié)點(diǎn)編號(hào)、網(wǎng)絡(luò)劃分以及邊界條件輸入等工作；在分析主程序模塊中，求解整體剛度方程；在后處理模塊中，實(shí)現(xiàn)結(jié)果顯示、數(shù)據(jù)輸出等工作。

對(duì)應(yīng)的有限元法的基本步驟：

（1）幾何域離散，獲得標(biāo)準(zhǔn)化的單元； （2）通過(guò)能量原理（虛功原理或最小勢(shì)能原理，獲得單元?jiǎng)偠确匠蹋?（3）單元的集成（裝配）； （4）處理位移邊界條件； （5）計(jì)算位移場(chǎng)；

（6）計(jì)算單元的其他物理量（應(yīng)力應(yīng)變）。

這幾步中，最核心的內(nèi)容是單元研究，具體包括：

（1）節(jié)點(diǎn)描述（不同坐標(biāo)系節(jié)點(diǎn)坐標(biāo)的變化）； （2）場(chǎng)描述（位移場(chǎng)，應(yīng)變場(chǎng)，應(yīng)力場(chǎng)，形函數(shù)）；

（3）單元?jiǎng)偠确匠蹋ɑ谀芰吭硗茖?dǎo)）。

需要說(shuō)明的是后文的四面體單元有限元方程的推導(dǎo)過(guò)程是基于等參單元的基本理論從局部坐標(biāo)（自然坐標(biāo)、體積坐標(biāo)）出發(fā)來(lái)推導(dǎo)四面體單元的剛度矩陣，因?yàn)檫@樣做比較規(guī)范自然，推導(dǎo)過(guò)程也適用于其他類型單元。但是因?yàn)樗拿骟w單元相對(duì)簡(jiǎn)單也可以直接從直角坐標(biāo)（全局坐標(biāo)）進(jìn)行推導(dǎo)，具體推導(dǎo)過(guò)程可參考清華大學(xué)曾攀老師的課程，直接從直角坐標(biāo)（全局坐標(biāo)）進(jìn)行推導(dǎo)的過(guò)程省去了等參單元雅各比矩陣呀等坐標(biāo)系映射的各種概念，理解起來(lái)相對(duì)容易。

公式（1）-（3）為彈性力學(xué)中三維空間彈性問(wèn)題的完整描述，分別是空間問(wèn)題的平衡方程、幾何方程、物理方程，這是我們推導(dǎo)有限元方程的基礎(chǔ)。這些公式會(huì)在后面的有限元方程推導(dǎo)過(guò)程中用到。

四面體單元的坐標(biāo)描述涉及了等參單元的概念，在有限元方法中，若要離散邊界為曲線或曲面的求解域，需要建立將形狀規(guī)則的單元變換為邊界為曲線或曲面的單元的方法，在有限元法中對(duì)應(yīng)此問(wèn)題所采用的變換方法是等參變換，即單元幾何形狀的變換和單元內(nèi)場(chǎng)函數(shù)采用相同數(shù)目的節(jié)點(diǎn)及相同的插值函數(shù)進(jìn)行變換。四面體單元的參數(shù)坐標(biāo)就是體積坐標(biāo)，體積坐標(biāo)的定義如圖2所示。

圖2 四面體單元體積坐標(biāo)的定義

參數(shù)坐標(biāo)系下的形函數(shù)等于對(duì)應(yīng)的體積坐標(biāo)，四個(gè)節(jié)點(diǎn)對(duì)應(yīng)四個(gè)形函數(shù)，如下式，

這樣就實(shí)現(xiàn)了自然坐標(biāo)系和物理坐標(biāo)系下的坐標(biāo)映射，如下式，

同樣形函數(shù)也可以用于物理場(chǎng)的插值，公式（6）是位移場(chǎng)的插值表達(dá)式。

有了位移場(chǎng)插值的表達(dá)式就可以通過(guò)公式（2）-（3）的幾何方程和物理方程推導(dǎo)應(yīng)變場(chǎng)和應(yīng)力場(chǎng)的有限元表達(dá)式，但是公式（2）的幾何方程中存在對(duì)于xy也就是物理坐標(biāo)系下的偏導(dǎo)數(shù)，代入公式（6）后也就是形函數(shù)對(duì)物理坐標(biāo)的偏導(dǎo)數(shù)，因?yàn)槲覀兊男魏瘮?shù)是在自然坐標(biāo)系下定義的，是關(guān)于ξ, η, ζ的函數(shù)，想要知道他對(duì)x，y的偏導(dǎo)，這里利用了鏈?zhǔn)角髮?dǎo)法則，即可建立如下式（7）所示的形函數(shù)對(duì)自然坐標(biāo)的物理坐標(biāo)偏導(dǎo)數(shù)的映射關(guān)系。中間的這個(gè)矩陣就叫Jacob矩陣。

如何求Jacob矩陣呢？將公式（5）代入公式（7）可以得到如下式（8）所示的自然坐標(biāo)偏導(dǎo)矩陣乘以一個(gè)坐標(biāo)矩陣，

得到j(luò)acob矩陣后求逆，再乘以自然坐標(biāo)偏導(dǎo)矩陣就可以得到形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)。如下式（9），

上述形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)的求解過(guò)程對(duì)應(yīng)的matlab代碼如下：

function [NDxyz,JacobiDET] = ShapeFunction(ElementNodeCoordinate)
%計(jì)算形函數(shù)及形函數(shù)對(duì)局部坐標(biāo)ksi eta zeta的導(dǎo)數(shù)
NDL = [-1 1 0 0;-1 0 1 0;-1 0 0 1];%3*4  [N1Dksi N2Dksi N3Dksi N4Dksi;N1Deta N2Deta N3Deta N4Deta……]
Jacobi = NDL*ElementNodeCoordinate;%計(jì)算雅可比矩陣3*4 4*3
JacobiDET = det(Jacobi);%計(jì)算雅可比行列式3*3 [DxDksi DyDksi DzDksi;DxDeta……
JacobiINV=inv(Jacobi);%對(duì)雅可比行列式求逆3*3
NDxyz=JacobiINV*NDL;%利用雅可比行列式的逆計(jì)算形函數(shù)對(duì)結(jié)構(gòu)坐標(biāo)的導(dǎo)數(shù)[DN1Dx DN2Dx DN3Dx;DN1Dy DN2Dy DN3Dy;……]
end

這樣求出形函數(shù)對(duì)物理坐標(biāo)的導(dǎo)數(shù)后就可以代入公式（2）幾何方程求出應(yīng)變場(chǎng)矩陣B，經(jīng)過(guò)能量原理的推導(dǎo)可以得到單元?jiǎng)偠染仃嚨谋磉_(dá)式，注意剛度矩陣的表達(dá)式是一個(gè)積分的運(yùn)算，由于被積函數(shù)較為復(fù)雜，如果代數(shù)積分進(jìn)行計(jì)算要消耗大量的計(jì)算資源，因此有限元理論中引入數(shù)值積分，即我們熟悉的高斯積分和hammer積分，對(duì)于直角坐標(biāo)系我們采用高斯積分，對(duì)于面積或者體積坐標(biāo)系我們采用hammer數(shù)值積分，具體這兩類積分的原理大家可以參考相關(guān)數(shù)值積分教材即可，這里我們只利用hammer數(shù)值積分的結(jié)論。四面體單元具體的數(shù)值積分公式如下：

其中，對(duì)應(yīng)的B矩陣如下式所示：

單元?jiǎng)偠染仃嚧_定好之后，將其按照自由度索引組裝到全局剛度矩陣中，組裝的核心思想就是確定好每個(gè)單元中的每個(gè)節(jié)點(diǎn)中的每個(gè)自由度在整體剛度矩陣中的位置即可。

整體剛度矩陣確定之后，就是荷載向量p的定義，這個(gè)很簡(jiǎn)單，只要依次確定每個(gè)節(jié)點(diǎn)每個(gè)自由度的外力荷載即可。

剛度矩陣和荷載向量確定后接下來(lái)就是邊界條件的施加，大家最熟悉的可能就是直接劃行劃列的方法，這個(gè)方法的弊端在于劃去自由度為零的行列之后，整體單元?jiǎng)偠染仃嚨木幪?hào)以及自由度編號(hào)都會(huì)改變，對(duì)于大規(guī)模計(jì)算非常不方便，影響計(jì)算效率。為了解決這個(gè)問(wèn)題出現(xiàn)了乘大數(shù)法和置一法，置一法只能解決零邊界問(wèn)題，乘大數(shù)法可以解決零邊界和非零邊界，應(yīng)用非常廣泛。邊界條件的施加這里我們采用乘大數(shù)法。針對(duì)邊界條件采用乘大數(shù)法施加邊界的方式為：在自由度上乘以一個(gè)很大的數(shù)a，然后在相應(yīng)的外力向量對(duì)應(yīng)的自由度上乘以，具體如下式（22）所示

為了驗(yàn)證乘大數(shù)法的可靠性，將對(duì)應(yīng)行列寫(xiě)成代數(shù)表達(dá)式的形式，如下式所示

考慮到遠(yuǎn)大于，所以公式（13）可以寫(xiě)成，

即可得到。可見(jiàn)乘大數(shù)法與實(shí)際邊界條件是等價(jià)的。

施加了邊界條件的剛度矩陣K與荷載矩陣p之后，可直接利用matlab中的高斯消去法 計(jì)算符“”，完成位移的求解，即u=Kp。

求出節(jié)點(diǎn)位移之后通常我們做受力分析時(shí)還會(huì)需要知道應(yīng)力應(yīng)變的分布情況，具體通過(guò)下述公式求解：

對(duì)應(yīng)的提取應(yīng)力和應(yīng)變的后處理代碼如下：

% 計(jì)算形函數(shù)導(dǎo)數(shù)
    [NDxyz, ~] = ShapeFunction(ElementNodeCoordinate);%[DN1Dx DN2Dx DN3Dx;DN1Dy DN2Dy DN3Dy;……]
    ElementNodeDisplacement=U(ElementNodeDOF);%12*1 節(jié)點(diǎn)位移列陣
    ElementNodeDisplacement=reshape(ElementNodeDisplacement,Dof,4);%3*4
    % 計(jì)算單元應(yīng)變 Strain3_3 3*3的應(yīng)變矩陣
    Strain3_3=ElementNodeDisplacement*NDxyz';%高斯積分點(diǎn)處應(yīng)變 3*4  4*3
    %把單元應(yīng)變矩陣改寫(xiě)成6*1
    Strain=[Strain3_3(1,1) Strain3_3(2,2) Strain3_3(3,3) ...
    Strain3_3(1,2)+Strain3_3(2,1)....
    Strain3_3(2,3)+Strain3_3(3,2) Strain3_3(1,3)+Strain3_3(3,1)]';
Stress(1:6,1) = D*Strain;%高斯積分點(diǎn)處應(yīng)變

利用繪圖命令Patch我們可以得到如圖1所示的應(yīng)力云圖。另外計(jì)算得到的變形和位移云圖如圖3-4所示。



圖3 變形前后的網(wǎng)格對(duì)比

圖4 位移云圖

審核編輯：劉清