譚 宸

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

0 引言

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

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

1 彈性問(wèn)題的有限單元格式

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

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

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

(1)

(2)

用笛卡爾張量符號(hào)可以進(jìn)一步將彈性力學(xué)平衡方程改寫(xiě)為:

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

(3)

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

(4)

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

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

(5)

改寫(xiě)成矩陣形式可得:

(6)

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

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

(7)

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

{ε}=ε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)

局部單元?jiǎng)偠染仃嘖:

(10)

代入B矩陣,進(jìn)而可以得到:

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

其中,

。

2 程序架構(gòu)

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

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

2.1 前處理程序

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

前處理程序的首要是生成節(jié)點(diǎn)坐標(biāo),包括節(jié)點(diǎn)的編號(hào)及相應(yīng)的點(diǎn)的坐標(biāo)值。節(jié)點(diǎn)位置信息定義為函數(shù)square_node_array,生成節(jié)點(diǎn)信息后,需要獲取每個(gè)單元的節(jié)點(diǎn)信息。節(jié)點(diǎn)的空間信息存儲(chǔ)在element數(shù)組中。為便于后續(xù)邊界條件的處理,方便使用,按照單元中節(jié)點(diǎn)逆時(shí)針編號(hào)進(jìn)行存儲(chǔ)記錄,在完成節(jié)點(diǎn)的計(jì)算、單元節(jié)點(diǎn)信息的輸入后,編寫(xiě)supportcond函數(shù)引入邊界條件。對(duì)于懸臂梁在端部受集中荷載而言,其荷載位移邊界條件的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 有限元計(jì)算分析程序

其中,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 后處理可視化程序

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

3 數(shù)值算例計(jì)算與驗(yàn)證

3.1 算例計(jì)算

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

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

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

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

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

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

3.2 方法驗(yàn)證

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

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

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

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

4 結(jié)論

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