結(jié)點(diǎn)矩形單元的有限元公式推導(dǎo)

1、六結(jié)點(diǎn)矩形單元的有限元公式推導(dǎo)一. d型6結(jié)點(diǎn)矩形單元2aX562b1.1單元位移模式及插值函數(shù)的構(gòu)造首先引入局部坐標(biāo)纟、其中 = - (一 151)耳=a b利用教材中Serendipity四邊形單元的構(gòu)造方法，假設(shè)開始只有4的個角結(jié)點(diǎn)，對應(yīng)這些結(jié)點(diǎn)的插值函數(shù)可以用雙一次拉格朗日多項(xiàng)式構(gòu)造，即凡=扛1 + )(1 + )N2 =1(1+0(1-77)N3=l(l-Xl-7)N4=1(1-)(1+7Z)當(dāng)增加結(jié)點(diǎn)5、結(jié)點(diǎn)6時，與它對應(yīng)的插值函數(shù)可表示為y方向2次，x方向1次拉格 朗日多項(xiàng)式的乘積，即他=*(1 + 歹)(1-2)他#(1-歹)(1- 2)此時77,(/= 1,2,3,4)在結(jié)點(diǎn)

2、5、結(jié)點(diǎn)6處不為0,為了使其繼續(xù)滿足這個要求，將英修正為A1M =心_號弘=_*（1 + 勺（一2 + 乃+ ）八1.I!N產(chǎn)N,7N嚴(yán)-屮-綁-於+亍+ p八1all4 =4-16 =-(1-)(-72-/7 + -)繼而可將位移可以證明這樣構(gòu)造的插值函數(shù)滿足叫=5.和為M = 1這些基本要求,1-1函數(shù)表示成結(jié)點(diǎn)位移的函數(shù)u = Ng + N2u2 + Nyq + N4u4 + N5u5 + Nj%v =+ A2v2 + Ny、+ N + N5* + A6v6矩陣形式為N、00N=Nae12應(yīng)變矩陣的推導(dǎo)/ 、exbxs0 / 、 uydyf xy )d d、內(nèi) 6x）d0=Lu = LN

3、ae = BaB = (BB.B,幾）oNtdxdaydxoNt01.3應(yīng)力矩陣的推導(dǎo)/o 1V0z、Xby= Ds =EV10X1 一芮J001-vI2丿S=DB= DBa=Sae1.4單元剛度矩陣的推導(dǎo)利用最小位能原理建立有限元方程，求得單元剛度矩陣為KeDBtdxdy其中BrDB =B；D(B B2)=oNtb：db*dx6oNiB；DBBDB,BDBB；DB2B；DB：BDB2b；db3 BDB. B；DB,001-vdN：dxoNtdNi+2VL dx dyJv 6叫込dNt dNj2dxdy dy(bl ia r1每 a+2叫+vr1一-i-ba D +12Do嚴(yán)+2 dx dx

4、v6N.級dx dxdN. 6N tQ2-v bA；2a)2 dy dx v6N, dNj其中A；=dNi Nj 蝕M dN ON：“ dN ON.B產(chǎn)JJ抽盤0% drfr dN- dNtD汀JJ d抽 -&l期的形函數(shù)的偏導(dǎo)數(shù)如下表BN、dN3W=_如_1)(+)典Y 1)(“dr4線=(+1)(+orj4蟲=+1)(-or4怨十1)(+3d/4帖 2二. b型6結(jié)點(diǎn)矩形單元12a2b2. 1單元位移模式及插值函數(shù)的構(gòu)造同樣引入局部坐標(biāo)纟、其中 = - (_1“G) /? = 7 (-171)a b假設(shè)開始只有4個角結(jié)點(diǎn),對應(yīng)這些結(jié)點(diǎn)的插值函數(shù)可以用雙一次拉格朗日多項(xiàng)式構(gòu)造, 即U=4(i

5、+g)(i+)N2 =1(1+0(1-77)N3=l(l-)(l-77)當(dāng)增加結(jié)點(diǎn)5時，與它對應(yīng)的插值函數(shù)可表示為y方向2次，x方向1次拉格朗日多項(xiàng) 式的乘積，當(dāng)增加結(jié)點(diǎn)6時，與它對應(yīng)的插值函數(shù)可表示為x方向2次，y方向1次拉格朗 日多項(xiàng)式的乘積，即他=*（1 + 歹）（1-九=4（1一鬥（1 + ）此時庖（：= 1,2,3,4）在結(jié)點(diǎn)5、結(jié)點(diǎn)6處不為0,為了使其繼續(xù)滿足這個要求，將英修 正為八1=-抽5-則6 蔦（ + 1）（ + 1）（ + -1）川2=八24弘=扌 + 1）（于_）八1”3=川3=才（1）N廠心-則6=?孑-刈+1）可以證明這樣構(gòu)造的插值函數(shù)滿足N廠和這些基本要求，繼而可

6、將位移1-1函數(shù)表示成結(jié)點(diǎn)位移的函數(shù)u = Ng + Njs + N3M3 + NjU + N、ig + Njsv = N“ + N2v2 + N3V3 + Nm + N5V5 + N6v6矩陣形式為0 他 0 N3 0 y 0 他 0 M 丿2. 2應(yīng)變矩陣的推導(dǎo)2. 3應(yīng)力矩陣的推導(dǎo)英中btdb =d0/ 、oxbx0d/ 、 uydy/巧丿dd6= Lu = LNae =BaB = (BB = LN,=B,adxdoyoNidxdNj6oNt彷( (TI V0/ XEXb,/)V10r1-/.rr001-vYW2丿&丿S=DB= DBa=Sae2. 4單元剛度矩陣的推導(dǎo)利用最小位能原理建

7、立有限元方程，求得單元剛度矩陣為KeDBtdxdyB；B；B,B；DB,Bj2DBB；DBb；db2B；DB：B；DB?b；db3bdb3B；DB、伽i0型、/1 V0dx0BjDB廠已dx勿V 100dN.1J 1-v20dNtoNt1-v労Idy&丿0 02dNtoNfl勿dx )其中+dx dx 2 dy dy67V dN _v6Nj 6Nv- +vdx dy 2 dy dx2 dx6勿dyb A 1-v A. +a !J2D b ,J% +VC Jib掃+5+ 2Du2 dx dx %2 i;-v b A TaM州=dNi込rr oN ON.Gj= JJ ddri.也77fr dN.

8、dN.Dj= JJ -、蝕其中形函數(shù)的偏導(dǎo)數(shù)如下表譽(yù)冷（2屮）（心）AN 1詬蔦（2 土 ）（屮3（宀）4 / 7心=扣 + 1）（2-1）d/4亠冶1）drf 4藥=扣7di 4空詁（1-孑）8 2三. a型6結(jié)點(diǎn)矩形單元與4結(jié)點(diǎn)矩形單元的計算比較3.1計算示例與參數(shù)為了驗(yàn)證上文中6結(jié)點(diǎn)矩形單元有限元公式的正確性，此處對圖示結(jié)構(gòu)用6結(jié)點(diǎn)矩形單 元與4結(jié)點(diǎn)矩形單元分別進(jìn)行受力分析。用6結(jié)點(diǎn)矩形單元進(jìn)行有限元分析的程序用C語言 自行編寫，而用4結(jié)點(diǎn)矩形單元進(jìn)行有限元分析則在ANSYS中進(jìn)行。圖示結(jié)構(gòu)為一正方形薄板，邊長為10m,左下角、右下角受到固泄約朿，左上角受到1N 的拉力，右上角收到1N的

9、壓力，共劃分為100個單元。彈性模量E.薄板厚度t均取為b泊松比取為0. 25o具體圖示如下。1ANF21314IS1G171919202112F3 27 201314:51:4411495fi6776IS34103112121231048576675M931021111202452566671)392101110H926BQ6bb64739291.100LD9118Z67066372919099LDB117276的62?1)00990107H6295435261?O79H83710611S29102516069?aB736L05113034150536B11R6951041133LX X90

10、3938373636343332N4結(jié)點(diǎn)矩形單元受力分析示意圖（利用ANSYS） 1A 1AAA 1 - ii J A iA jia AAA ,A 1 i ii 1 i AA JA 1a型6結(jié)點(diǎn)矩形單元受力分析示意圖32計算結(jié)果及對比利用ANSYS軟件訃算的位移結(jié)果如下圖所示:DISPlACaeirTSIEP-1 5UB =1TIME=12MX =17.945FEB 27 20131S：44:24由于使用自己編寫的程序計算的位移結(jié)果數(shù)據(jù)量太大,所以在此僅列出若結(jié)點(diǎn)的位移 值并與ANSYS的計算結(jié)果相比較。典型節(jié)點(diǎn)如下圖所示。在典型結(jié)點(diǎn)處，使用6結(jié)點(diǎn)單元和4結(jié)點(diǎn)單元求得的位移結(jié)果的比較見下表。結(jié)

11、點(diǎn)編號6結(jié)點(diǎn)X方向位移4結(jié)點(diǎn)(ANSYS)誤差Y方向位移4結(jié)點(diǎn)誤差113. 87112. 6891. 18213. 87112. 6891. 18227. 49186. 72980. 7624. 91454. 56150. 35334.20213. 69630. 50581.44311. 32400. 119141.42791. 15840. 2695-1.3311-1.2139-0.11725-0.83084-0.883780. 05294-3. 4081-3.0521-0. 356613. 87112. 6891. 182-13. 871-12. 689-1. 18277. 49186.

12、72980. 762-4.9145-4.5615-0. 35384.20213. 69630. 5058-1.4431-1. 3240-0.119191.42791. 15840. 26951. 33111.21390. 117210-0.83084-0.883780. 052943. 40823. 05210. 35613. 3結(jié)果分析從計算得到的位移結(jié)果可以發(fā)現(xiàn)，使用6結(jié)點(diǎn)與4結(jié)點(diǎn)單元進(jìn)行分析得到的位移誤差較 小，這驗(yàn)證了本文之前所推導(dǎo)得到的有限元公式的正確性。但是，兩者之前仍然具有一泄的 誤差。這可能是由于劃分單元的個數(shù)過小造成的。四、6結(jié)點(diǎn)單元C語言計算程序#include#inclu

13、deinclude FILE *fpElist, *fpNlist, *fpDlist, *fpFlist,*fpOutput;FILE *fpNodeData, *fpElemData;/輸出文件，用于 ANSYS 的輸入define NMAX 1000/最大單元數(shù)減一double t=l;/薄板厚度double Ey=l;/楊氏模量double v=0. 25;/泊松比int i,j,k;/循環(huán)變量double dbTemp, dbTempl;/臨時變量int NumOfNode;/ 點(diǎn)總數(shù)int ENMAX+1 7, NumOfElem;/單元信息矩陣，單元總數(shù)double KCNMAX

14、+1 ENMAX+1, kk13 13;總體剛度矩陣，單元剛度矩陣int嘰NMAX*3;/節(jié)點(diǎn)位移約束矩陣，1表示位移為0, 0表示位移不受約朿；排列順 序?yàn)椋汗?jié)點(diǎn)1的X方向，節(jié)點(diǎn)1的Y方向，節(jié)點(diǎn)2的X方向。double FNMAX*3;/節(jié)點(diǎn)荷載矩陣，排列順序?yàn)椋汗?jié)點(diǎn)1的X方向受力，節(jié)點(diǎn)1的Y 方向受力，節(jié)點(diǎn)2的X方向受力“。int Num0fVari=0,VENMAX*3;/自由變量總數(shù)，自由變量-所有變量double UNMAX*3;/節(jié)點(diǎn)位移變崑 結(jié)果double aa_66=1. /5,1. /30,-1./30,-1./5,-1./5,1./5, 1. /30,1./5,-1./5

15、,-1./30,-1./5,1./5, -1. /30,-1./5,1./5,1./30,1./5,-1./5, -1. /5,-1./30,1./30,1./5,1./5,-1./5, -1. /5,-1./5,1./5,1./5,8./15,-8./15, 1. /,1./5,-1./5,-1./5,-& /15,8./15;double bb66=1. /4,5. /12,5. /12,1./4,-1./3,-1./3, -5. /12,-1. /4,-1./4,-5. /12,1./3,1./35. /12,1./4,1/4,5. /12,-1./3,-1. /3-1. /4,-5./

16、12,-5.12,-1./4,1./3,1./31./3,-1./3,-1/3,1./3,0,0,-1./3,1./3,1/3,-1./3,0,0;double cc.66 = 1./4,-5./12,5./12,-1./4,1./3,-1. /35. /12,-1./4,1./4,-5./12,-1./3,1./35. /12,-1./4,1./4,-5./12,-1./3,1./31./4,-5./12,5./12,-1./4,1./3,-1. /3-1. /3,1./3,-1./3,1./3,0,0,-1. /3,1./3,-1./3,1./3,0,0；19. /9,13. /9,13.

17、/I&19./I&-16. /9,13. /9,19. /9,19./I&13./I&-16. /9,13. /1&19. /18,19. /9,13. /9,-8 /9,19. /1&13./I&13. /9,19. /9,-& /9,-16. /9,-16. /9,-& /9,-8. /9,16. /9,-8. /9,一& /9,-16. /9,-16. /9,& /9,double dd66=-8. /9,/9, -16./9, -16. /9,8. /9,16. /9;double a_b=l;/矩形單元的長寬比a/bvoid mainO/建立單元信息矩陣，100個單元，11*21個結(jié)

18、點(diǎn)Num0fElem=100;Num0fNode=U*21;for(i=0;i=9;i+)for(j=l;j=10;j+)Ei*10+jl=22*i+j+23;Ei*10+j2=22*i+j+l;Ei*10+j3=22*i+j；Ei*10+j4=22*i+j+22;Ei*10+j5=22*i+j+12;Ei*10+j 6=22*i+j+U;for (i=l; i=NumOfElem; i+)printfC%d) %d %d %d %d %d %dn, i, Ei 1, Ei 2, Ei 3, Ei 4, Ei 5, Ei6);計算單元剛度矩陣 for(i=l;i=6;i+)for(j=l;j=

19、6;j+)kk2*i-l 2*j-l=aai-l j-l+(l-v)/2*ddi-l j-1;kk2*i-l 2*j=v*bbi-l j-l+(l-v)/2*cci-l j-1;kk2*i2*j-l=v*cci-lj-l+(l-v)/2*bbi-lj-1; kk2*i2*j=ddi-lj-l+(l-v)/2*aai-lj-1; printfC 單剛 rT);for(i=l;i=12;i+)for(j=l;j=12;j+)kki j=kki j*Ey/(l-v*v); printf (z，%lf 、kki j);printf (z/nz，);輸出單元剛度矩陣的結(jié)果文件if (fpOutput=f

20、openCDanGang txt, w+)=NULL) printf(Cannot open this filen); exit(O);for(i=l;i=12;i+)for(j=l;j=12;j+)fprintf(fpOutput, %741f ”、kkij); fprintf(fpOutput, n);fclose(fpOutput);/根據(jù)單元剛度矩陣裝配整體剛度矩陣for(k=l;k=NumOfElem;k+)for(i=l;i=6;i+)for(j=l;j=6;j+)KEki*2-lEkj*2-l+=kki*2-lj*2-l;KEk i*2-lEk j*2 +=kki*2-lj*2

21、;KEk i*2 Ek j*2-l+=kki*2 j*2-l;KEki*2 Ekj*2 +=kki*2 j*2 ;輸岀整體剛度矩陣的結(jié)果文件if (fpOutput=fopen(ZhengGang txt, w+)=NULL)printf (Cannot open this file n);exit(O);for (i=l; i=2*Num0fNode; i+)for(j=l;j=2*Num0fNode;j+)fprintf (fpOutput, /z%7. 41f ”、Ki j);fprintf(fpOutput, n);fclose(fpOutput);設(shè)宜盯點(diǎn)位移約束矩陣D 1=1;D 2二 1; D221*2-1=1;D221*2=1;/設(shè)巻肖點(diǎn)荷載矩陣Fll*2-l=l;F231*2-1二-1;/*/解線性方程