基于MATLAB平臺的懸臂梁靜力彈性分析

譚 宸

(同濟大學(xué)土木工程學(xué)院,上海 200092)

0 引言

有限單元法是當今工程分析中獲得廣泛應(yīng)用的數(shù)值方法[1],采用有限元方法可以求解具有復(fù)雜的幾何邊界的彈性問題。曹子龍等[2]建立了求解圓孔薄板彈性問題的雜交應(yīng)力單元,推導(dǎo)出圓孔薄板問題的極坐標下應(yīng)力插值矩陣,可較高精度的求解孔邊附近應(yīng)力。陳立勝等[3]采用有限單元法初步對彈性平面問題進行了初步嘗試。張明哲等[4]開發(fā)的計算軟件平臺有助于學(xué)生很好的理解在彈塑性問題中改進無單元Galerkin方法的應(yīng)用。目前,伴隨著有限元軟件功能的不斷完善,對于剛接觸有限元軟件的學(xué)生而言,往往對其一知半解,在有限元分析中著重于物理參數(shù)的調(diào)節(jié)[5]。

為了加深學(xué)生們對于有限元方法基本原理的理解,本文通過基于MATLAB平臺,編寫結(jié)構(gòu)化的有限元代碼,可以應(yīng)用于基本的結(jié)構(gòu)工程問題。該代碼的結(jié)構(gòu)易于理解,并可由學(xué)生輕松擴展。本文開發(fā)的代碼計算懸臂梁變形結(jié)果與商業(yè)有限元軟件ABAQUS計算結(jié)果一致,開發(fā)的一個結(jié)構(gòu)良好、易于使用的計算懸臂梁結(jié)構(gòu)彈性變形的有限單元分析程序具有很好的適用性。

1 彈性問題的有限單元格式

工程中許多問題,通常是以未知場函數(shù)在域內(nèi)應(yīng)滿足的偏微分方程形式提出。根據(jù)彈性力學(xué)基本方程和與之等效的變分原理,可以從基本方程的強形式推導(dǎo)出有限元格式,進而可獲得局部單元剛度矩陣。

本小節(jié)以二維平面域中彈性問題為例,考慮邊界條件,對其偏微分方程的有限元格式進行推導(dǎo)。從面積域中取一無限小單元,其受力示意圖如圖1所示。

從圖1可以看出,該單元的平衡方程可表示為:

(1)

(2)

用笛卡爾張量符號可以進一步將彈性力學(xué)平衡方程改寫為:

σji,j+bi=0 (i,j=1,2)

(3)

通過加權(quán)余量法可寫出平衡微分方程的等效積分弱形式。等效積分弱形式可以通過分部積分得到:

(4)

作為平衡微分方程的等效積分弱形式,在導(dǎo)出過程中并未涉及物理方程,因此,不僅僅可以適用于線彈性問題,還可以用于非線性彈性或者彈塑性問題。

將權(quán)函數(shù)wi替換成δui,可以得到有限單元法的控制方程:

(5)

改寫成矩陣形式可得:

(6)

本文采用4結(jié)點矩形單元來處理二維平面彈性力學(xué)問題,4結(jié)點矩形單元如圖2所示。

對于每個角點,構(gòu)造出它的形函數(shù)為:

(7)

進一步,根據(jù)彈性力學(xué)中位移與應(yīng)變的關(guān)系,可將單元應(yīng)變用結(jié)點位移向量來進行表示:

{ε}=εxεyγxyì?í????üty????=N1(x,y)x0N2(x,y)x0N3(x,y)x0N4(x,y)x00N1(x,y)y0N2(x,y)y0N3(x,y)y0N4(x,y)yN1(x,y)yN1(x,y)xN2(x,y)yN2(x,y)xN3(x,y)yN3(x,y)xN4(x,y)yN4(x,y)xé?êêêêêêêù?úúúúúúúüty??????????????????????????????????[B] u1 v1 u2 v2 u3 v3 u4 v4u1v1u2v2u3v3u4v4ì?í???????????üty???????????(8)

將形函數(shù)代入可得梯度矩陣[B]:

(9)

局部單元剛度矩陣K:

(10)

代入B矩陣,進而可以得到:

[K]= u1 v1 u2 v2 u3 v3 u4 v4k11k12k13k14k15k16k17k18k21k22k23k24k25k26k27k28k31k32k33k34k35k36k37k38k41k42k43k44k45k46k47k48k51k52k53k54k55k56k57k58k61k62k63k64k65k66k67k68k71k72k73k74k75k76k77k78k81k82k83k84k85k86k87k88é?êêêêêêêêêêêù?úúúúúúúúúúúu1v1u2v2u3v3u4v4(11)

其中,

。

2 程序架構(gòu)

該程序基于MATLAB編程語言開發(fā),依托于有限單元法解決彈性問題,可以高效便捷的對懸臂梁受集中荷載作用下的應(yīng)力應(yīng)變進行分析。程序通過先對幾何域進行單元網(wǎng)格劃分,建立節(jié)點和對單元進行編號,隨后輸入邊界條件等模型數(shù)據(jù),對其進行采用有限單元法的彈性分析后得到處理結(jié)果,最后將結(jié)果云圖以可視化的界面進行展示。程序可以分為前處理程序、有限元彈性分析計算程序和可視化后處理程序三個部分,其流程圖和MATLAB函數(shù)如圖3所示。

前處理程序主要是通過建立幾何模型,生成節(jié)點和網(wǎng)格,將分析域離散化,組建全局坐標系統(tǒng),同時引入邊界條件。輸入?yún)?shù)后預(yù)處理必備的結(jié)構(gòu)參數(shù)值,前處理關(guān)鍵環(huán)節(jié)需要定義節(jié)點、生成每個單元坐標信息、輸入載荷信息和邊界條件,為有限元法分析做準備。計算分析程序進行有限單元法的計算,最主要的是進行矩陣計算和求解。計算分析程序通過求解全局剛度矩陣獲得每個節(jié)點的位移{u},進而評估每個單元上的應(yīng)力和應(yīng)變。后處理程序模塊對每個單元建立節(jié)點位移和應(yīng)力向量,進而繪制加載后模型的位移、應(yīng)力云圖,實現(xiàn)計算結(jié)果的可視化。

2.1 前處理程序

前處理程序是數(shù)值計算中非常重要的環(huán)節(jié),影響后續(xù)數(shù)值計算的精度和效率。前處理可分為如下部分:輸入材料形狀大小、材料性質(zhì)參數(shù)、生成節(jié)點坐標等單元信息以及載荷和邊界條件定義等。

前處理程序的首要是生成節(jié)點坐標,包括節(jié)點的編號及相應(yīng)的點的坐標值。節(jié)點位置信息定義為函數(shù)square_node_array,生成節(jié)點信息后,需要獲取每個單元的節(jié)點信息。節(jié)點的空間信息存儲在element數(shù)組中。為便于后續(xù)邊界條件的處理,方便使用,按照單元中節(jié)點逆時針編號進行存儲記錄,在完成節(jié)點的計算、單元節(jié)點信息的輸入后,編寫supportcond函數(shù)引入邊界條件。對于懸臂梁在端部受集中荷載而言,其荷載位移邊界條件的matlab代碼如下:

function [topEdge,topEdge1,dispNodes,dispNodes1,leftNodes1]=supportcond(numx,numy)

nnx=numx+1;

nny=numy+1;

uln=nnx*(nny-1)+1;

urn=nnx*nny;

lrn=nnx;

lln=1;

topEdge=[uln:1:(urn-1);(uln+1):1:urn]′;

topEdge1=topEdge;

botEdge=[lln:1:(lrn-1);(lln+1):1:lrn]′;

rightEdge=(lrn:nnx:(urn))′;

botNodes=unique(botEdge);

topNodes=unique(topEdge);

rightNodes=unique(rightEdge);

leftNodes=rightNodes-(nnx-1);

dispNodes=botNodes;

rightNodes1=rightNodes(2:end);

leftNodes1=leftNodes(2:end);

dispNodes1=leftNodes;

end

2.2 有限元計算分析程序

其中,B矩陣的matlab代碼如下:

function Bfem=Bmatrix(pt,iel)

global node element

sctr=element(iel,:);

nn=length(sctr);

[N,dNdxi]=shape_func(pt);

J0=node(sctr,:)′*dNdxi;

invJ0=inv(J0);

dNdx=dNdxi*invJ0;

Bfem=zeros(3,2*nn);

Bfem(1,1:2:2*nn)=dNdx(:,1)′;

Bfem(2,2:2:2*nn)=dNdx(:,2)′;

Bfem(3,1:2:2*nn)=dNdx(:,2)′;

Bfem(3,2:2:2*nn)=dNdx(:,1)′;

end

2.3 后處理可視化程序

后處理程序獲取有限元計算分析程序得到的位移等信息,并重新組織成標準后處理結(jié)果文件的形式,利用matlab可視化技術(shù)再現(xiàn)網(wǎng)格信息,將變形前后的網(wǎng)格展示出來。同時,分析計算結(jié)果,如位移云圖、應(yīng)力云圖的結(jié)果以圖像形式顯示出來。本系統(tǒng)的后處理程序包括數(shù)據(jù)處理程序、圖像顯示程序。

3 數(shù)值算例計算與驗證

3.1 算例計算

對二維受集中力的懸臂梁采用彈性問題的有限單元法進行小變形數(shù)值分析。通過對該懸臂梁進行自主編程建模,建模核心代碼流程如第2小節(jié)所示。利用前處理程序模塊在梁內(nèi)布置節(jié)點,進行單元編號,記錄每個單元的節(jié)點信息。將輸入節(jié)點信息和荷載位移邊界條件輸入到有限元分析計算程序中,進而計算獲得節(jié)點的位移和應(yīng)力、應(yīng)變等數(shù)據(jù)。將所得結(jié)果與ABAQUS軟件分析結(jié)果進行比較,驗證了本文提出的彈性問題有限單元法程序的有效性,可進一步加深學(xué)生們采用有限元法分析問題編制程序的理解。

懸臂梁自由端受集中荷載示意圖如圖4所示。梁的幾何尺寸為L=10 m,D=1 m,梁的厚度為t=0.1 m。梁受集中荷載P=0.1 N。材料的彈性模量為E=100 kPa,泊松比為v=0.25,按照平面應(yīng)力問題進行求解。

采用前處理程序進行結(jié)點生成,單元網(wǎng)格劃分,水平劃分單元數(shù)為100,豎向劃分單元數(shù)為10,結(jié)果如圖5所示。

將節(jié)點信息和荷載條件輸入到有限元分析計算程序中進行計算,獲得變形后的懸臂梁有限單元網(wǎng)格如圖6所示。

懸臂梁的位移云圖如圖7所示。

從圖7中可以看出,右端部頂處豎向結(jié)點位移為0.039 8 m。水平應(yīng)力σxx與豎向應(yīng)力σyy云圖見圖8。

3.2 方法驗證

基于ABAQUS有限元的數(shù)值模擬方法是檢驗本文程序編制的正確性的有效途徑之一。本小節(jié)通過建立與3.1算例相同的懸臂梁結(jié)構(gòu)有限元模型,并對模型賦予材料參數(shù)、施加集中荷載邊界條件,驗證本文方法的正確性。

在ABAQUS中,建立的懸臂梁三維結(jié)構(gòu)幾何參數(shù)為L=10 m,D=1 m,t=0.1 m。懸臂梁材料彈性模量E=100 kPa,泊松比為v=0.25。有限元模型如圖9所示。在有限元模型中,所施加的載荷為力載荷。模型網(wǎng)格劃分如圖10所示。

ABAQUS計算獲得的懸臂梁的位移云圖如圖11所示,從圖11中可以看出,右端部頂處豎向結(jié)點位移為0.040 m。水平應(yīng)力σxx與豎向應(yīng)力σyy云圖見圖12。

通過將本文編寫的懸臂梁有限單元法程序與ABAQUS計算的結(jié)果進行對比,可以看出,所開發(fā)的基于MATLAB平臺的懸臂梁彈性分析程序可以較好的預(yù)測梁的小變形,所得計算結(jié)果與有限元軟件計算結(jié)果一致。

4 結(jié)論

本文以求解懸臂梁受集中荷載作用下變形為算例,介紹了有限單元法求解彈性問題的算法。通過將商業(yè)有限元分析軟件結(jié)果與本文程序編制結(jié)果對比,兩者一致,驗證了本文程序的有效性。該程序基于MATLAB平臺編制而成,程序簡單實用,計算效率比采用有限元軟件ABAQUS分析彈性問題的效率更高。本文展示的有限元求解彈性問題的方法通用性強,通過該程序,可以進一步加深學(xué)生們采用有限元法分析問題編制程序的理解,具有重要意義和參考價值。通過對該程序的學(xué)習,可以舉一反三,在該程序的基礎(chǔ)上修改荷載和位移邊界條件,進一步求解其他的彈性問題。