雙線性四邊形等參單元有限元程序

1、I#|&門|杵if|參|azmedorz|匕1対150氐包也夠神9京怪比22J?；?0.13雙線性四邊形等參單元有限元程序本程序采用matlab編寫。程序加載由用戶提供的前處理數(shù)據(jù)，包括網(wǎng)格數(shù)據(jù)和載荷數(shù)據(jù)。采用直接的數(shù)值運(yùn)算和matlab符號運(yùn)算兩種方法（可選擇）生成單元剛度矩陣。自動集成結(jié)構(gòu)剛度矩陣，選用直接解法求解線性方程組，解出節(jié)點(diǎn)位移。后處理過程中，程序計算了節(jié)點(diǎn)應(yīng)力值，對于共享節(jié)點(diǎn)的應(yīng)力程序采用各個單元計算值的平均，程序同時給出了單元最佳應(yīng)力點(diǎn)的應(yīng)力值！1使用說明用戶需要在目錄中給定的文件中按照既定的格式給出必要的前處理數(shù)據(jù)（后面有詳細(xì)說明）。然后將matlab的工作目錄設(shè)為我們提供

2、的“雙線性四邊形等參單元程序”目錄，然后在命令行中輸入main。并按下回車鍵即可。如圖1所示：lorke*IVdue荷“a|Tm-ITftcWote=genXotoStxess(minXotelocWoLJsairiO0.025kzioeo0.3k-kK=aeeeatd8X(nijriXotejlocWote,rE,nFEclcclearEbX03-9-3上午10:01%notebocic屛-04A3上午10:11lcT-04A3上午10:1221_1多卻|aotMxrbr|gV/JLA&f汗珂V二門LXB|二1_3觀?！榜?|JF惱科l左|函2圖1程序運(yùn)行說明1.1輸入數(shù)據(jù)程序的輸入數(shù)據(jù)需要

3、在“雙線性四邊形等參單元程序”目錄中給定的文件中完成。每一個.txt文件對應(yīng)一個矩陣，所以數(shù)據(jù)的輸入必要嚴(yán)格按照matlab矩陣加載文件的格式完成，即文件的一行對應(yīng)于矩陣行，數(shù)據(jù)之間用空格或者逗號隔開，分號或者換行符表示進(jìn)入矩陣另一列的輸入。如圖2:所示輸入的矩陣陣為一4*4的矩陣:q245、235656S9E,vp,u)3主要函數(shù)簡介變量說明：locNote為節(jié)點(diǎn)坐標(biāo)，KE為單元剛度矩陣，K為總剛度矩陣，nuniNote為單元節(jié)點(diǎn)號逆時針，11E為單元數(shù)，nF為節(jié)點(diǎn)自由度數(shù)，E為彈性模量，v為泊松比，p為平面應(yīng)力和平面應(yīng)變的標(biāo)志，p=l為平面應(yīng)力，p=0為平面應(yīng)變，h待求結(jié)構(gòu)厚度。Ev,h,

4、p,locNotenuinNote,dB=dateInput0加載輸入數(shù)據(jù)，返回彈性模量E,泊松比v等。kK=assembleK(numNoteocNote,nE,nFEvhp)為總剛度矩陣集成包括單元剛度矩陣計算子函數(shù)eleiriK=geiiKE(numNote,locNoteriurnElementE,v,hp)和elemK=genKE2(numNote,locNote,nurnElement,E,v,h,p)。返回結(jié)構(gòu)總剛度矩陣。elemK=genKE(numNoteocNote,numElementEvhp)采用直接的數(shù)值計算的方法求解單元剛度矩陣，并返回。eleinK=geiiECE

5、2(nuiYfbTote,locNote,nui-nElement,E,v,h,p)采用符號運(yùn)算求解單元剛度矩陣，并返回。p=loadF(numNoteocNote,nF,h)計算總體載荷列陣，并返回。K,P=modityK(pP,kK,nF,dB)采用乘大數(shù)的方法引入位移邊界條件，消除剛度矩陣奇異性，返回修正后的剛度矩陣和載荷列陣。strNoteenNoteStress(numNoteocNote,nE，nFEvpu)計算節(jié)點(diǎn)應(yīng)力，采用圍繞該節(jié)點(diǎn)的單元計算應(yīng)力的平均值作為該節(jié)點(diǎn)的應(yīng)力值shCent=genCentShess(nuiriNotelocNote,nE,nFE,v,p,u)計算單元

6、最佳應(yīng)力點(diǎn)即局部坐標(biāo)(0,0)點(diǎn)的應(yīng)力值。4計算原理4.1形狀函數(shù)4節(jié)點(diǎn)4邊形等參元采用雙線性插值函數(shù)，對如圖6所示的母單元，插值函數(shù)為：圖64節(jié)點(diǎn)4邊形母單元4我們將位移表示為節(jié)點(diǎn)位移的插值形式，用矩陣寫為:心N(6)其中：。=址*1，處川2上3川3川4川47N=N00NN200N2N30N40N300N44i=i我們選用相同的插值函數(shù)做等參變換，即:(2)等參單元中插值函數(shù)是基于母單元的局部坐標(biāo)的，因此需要進(jìn)行局部坐標(biāo)下微分與整體坐標(biāo)微分的轉(zhuǎn)化，即：（吋、=Rdx2丿J為雅可比矩陣:八丄(X3-X4)(l+)+(X2-X1)(l-?7)&2呼)(1帀旳3呼4)(1+)Al丿124(x4-x

7、1)(i-)+(x3-x2)(i+)(y4-yi)(i)+(y3-y3)(i+)J21丿22將（2）式代入得:方程（3）逆過來寫為:(5)一站-期(6)(6)或者:dx1丿22-丿12=H-厶Jn生丿2丿4.2單元剛度矩陣及整體剛度矩陣的集成4節(jié)點(diǎn)4邊形單元剛度矩陣從單元應(yīng)變能而來:(8)h為結(jié)構(gòu)厚度。應(yīng)變矩陣可以表示為位移的函數(shù):/、dudx勿duOv(9)采用前面的微分轉(zhuǎn)化公式(7)將可表示為局部坐標(biāo)的微分形式:=(10)這里：厶進(jìn)一步，我們有_丿12000_丿21厶丿11J竝-丿12.(H)y7dz/一引-av-avYfI7-101-耳01+帀0-1-710-r01+歹0r40乃-100

8、1+刁0.010T歹01+孑00l(20-1-7rvsV1=Ga(12)將式(11)、(12)代入式(10)得:(13)e=AGa=Ba到此我們已將應(yīng)變表示為單元節(jié)點(diǎn)位移的函數(shù)。進(jìn)一步我們有應(yīng)力矩陣:(14)a=DBa這里D為彈性矩陣。將(13)和(14)代入(8)式得:TOC o 1-5 h z111Ue=fJBYDBjdra=-afKJa(15)2-1-12這里單元剛度矩陣Ke為：切卩|d?d(16)-1-1至此，單元剛度矩陣的計算已經(jīng)完成。函數(shù)eleniK=geiiKE(numNote,locNoteiiumElement,E,v,h,p)就是嚴(yán)格按照這個過程寫的的程序。至于函數(shù)elem

9、K=genKE2(nuniNoteJocNotejiumElement,E,v,h,p)直接采用符號運(yùn)算，用插值函數(shù)的倒數(shù)將B表示出來，代入(16)即可求出單元剛度矩陣。具體的過程以及E的最終表示式程序里面有明確的表示，這里不再重復(fù)。整體剛度矩陣的集成采用對號加入的方法進(jìn)行，這里沒有什么深入的理論，具體的操作程序中有詳細(xì)過程，這里不作詳細(xì)的說明。4.3載荷計算一般的工程計算中，重力等體力是不考慮的，因此我們這里沒有涉及體力的計算，事實(shí)上體力的計算相對比較簡單可以按照下式計算。11巧jJJmr門卩出切(17)-1-1我們對表面壓力進(jìn)行了計算，工程應(yīng)用中，主要涉及均勻載荷和線性漸變載荷，對于高次的

10、函數(shù)載荷涉及很少，我們這里就只處理均勻載荷和線性漸變載荷。(18)表面壓力按照如下公式進(jìn)行:，緲伽第+(訥切對于4節(jié)點(diǎn)4邊形單元由于各邊都是直線，()2+卡項實(shí)際上就是單元邊界的長度，可以直接用兩個端點(diǎn)的坐標(biāo)表示為：(19)(拿+(字)2丄1-朋+()_防卡drdr(20)T=對于線性變化的載荷我們只需要用戶輸入單元邊界兩個端點(diǎn)的壓力P?,p：(以X方向為例)即可代入計算，當(dāng)Px1=Px2時即為均勻載荷。由于(18)式中插值函數(shù)在邊界上都是線性的，因此我可以直接在單元邊界上建立局部坐標(biāo)，選取一維線性插值函數(shù)Nl=l-s/L,N2=s/L,這樣單元邊界表面壓力載荷矩陣T就可以表示為：Nlp；+N

11、2p；.將(19)、(20)式代入(18)并采用局部坐標(biāo)系下的插值函數(shù)，得到：N00NN20N2Nlp+N2p；Np+N2p：(21)(18)(18)其中1為單元邊界長度:/=(xl-x2)2+(yl-y2)2p對(21)式積分求解最后得:心中+護(hù)訓(xùn)+護(hù)111.111.評+尹；評+護(hù)(22)(18)到此一個單元邊界的等效載荷計算完成，按照(22)式一次計算所有的單元邊界并集成到總體等效載荷列陣中即可。我們的程序就是按照(22)式完成的。4.4線性方程組求解及后處理這里沒有采用先進(jìn)的方法求解線性方程組，直接采用matlab命令:a=KP進(jìn)行求解。后處理過程主要是求節(jié)點(diǎn)應(yīng)力以及最佳應(yīng)力點(diǎn)應(yīng)力值。節(jié)

12、點(diǎn)應(yīng)力的計算按照(14)式進(jìn)行，采用圍繞節(jié)點(diǎn)的所有單元計算的平均值作為節(jié)點(diǎn)應(yīng)力值！5算例作為一個簡單的演算，這里計算一單向均拉平板的應(yīng)力場。板長lm,寬0.25111,厚0.025m、長度方向收3MPa的均勻拉力。材料彈性模量210e6,泊松比0.3。取長度方向上一半來計算，長度方向劃分兩個單元，寬度方向一個單元,5612圖7算例示意圖如圖7所示。輸入文件見程序中默認(rèn)的設(shè)置!計算結(jié)果如圖8所示：圖82單元運(yùn)行結(jié)果可以發(fā)現(xiàn)，雖然節(jié)點(diǎn)應(yīng)力在2、4節(jié)點(diǎn)處偏離理論解較大，但是單元最佳應(yīng)力點(diǎn)的結(jié)果和理論解吻合的的非常好。進(jìn)一步采用4個單元，在高度上也劃分兩個單元，計算結(jié)果如圖9所示。圖94單元計算結(jié)果同樣應(yīng)力在最佳應(yīng)力點(diǎn)和理論結(jié)吻合的非常好，但是其它在節(jié)點(diǎn)2、4、5處偏離較大。因此可以看出4節(jié)點(diǎn)4邊形等參單元由于插值函數(shù)的為線性的，所以精度較差，實(shí)際的計算中需要較