1、精選優質文檔-傾情為你奉上八節點等參單元平面有限元分析程序土木工程學院2011.2專心-專注-專業目錄計算力學課程大作業1. 概述通常情況下的有限元分析過程是運用可視化分析軟件(如ANSYS、SAP等)進行前處理和后處理，而中間的計算部分一般采用自己編制的程序來運算。具有較強數值計算和處理能力的Fortran語言是傳統有限元計算的首選語言。隨著有限元技術的逐步成熟,它被應用在越來越復雜的問題處理中，但在實際應用中也暴露出一些問題。有時網格離散化的區域較大,而又限于研究精度的要求，使得劃分的網格數目極其龐大，結點數可多達數萬個，從而造成計算中要運算的數據量巨大，程序運行的時間較長的弊端，這就延長

2、了問題解決的時間，使得求解效率降低。因為運行周期長，不利于程序的調試，特別是對于要計算多種運行工況時的情況；同時大數據量處理對計算機的內存和CPU 提出了更高的要求，而在實際應用中，單靠計算機硬件水平的提高來解決問題的能力是有限的。因此，必須尋找新的編程語言。隨著有限元前后處理的不斷發展和完善，以及大型工程分析軟件對有限元接口的要求，有限元分析程序不應只滿足解題功能，它還應滿足軟件工程所要求的結構化程序設計條件，能夠對存儲進行動態分配，以充分利用計算機資源，它還應很容易地與其它軟件如CAD 的實體造型，優化設計等接口。現在可編寫工程應用軟件的計算機語言較多，其中C語言是一個較為優秀的語言，很容

3、易滿足現在有限元分析程序編程的要求。C語言最初是為操作系統、編譯器以及文字處理等編程而發明的。隨著不斷完善，它已應用到其它領域，包括工程應用軟件的編程。近年來，C語言已經成為計算機領域最普及的一個編程語言，幾乎世界上所有的計算機都裝有C的編譯器，從PC機到巨型機到超巨型的并行機，C與所有的硬件和操作系統聯系在一起。用C 編寫的程序，可移植性極好，幾乎不用作多少修改，就可在任何一臺裝有ANSI、C編譯器的計算機上運行。C既是高級語言，也是低級語言，也就是說，可用它作數值計算，也可用它對計算機存儲進行操作。2. 編程思想本程序采用C語言編程，編制平面四邊形八節點等參元程序，用以求解平面結構問題。程

4、序采用二維等帶寬存儲整體剛度矩陣，乘大數法引入約束，等帶寬高斯消去法求解位移。在有限元程序中，變量數據需賦值的可分為節點信息，單元信息，載荷信息等。對于一個節點來說，需以下信息：節點編號（整型），節點坐標（實型），節點已知位移（實型），節點載荷（實型），邊界條件（實型）和節點溫度（實型）等。同樣，對于一個單元來說，需以下信息：單元的節點聯接信息（整型），單元類型信息（桁架、梁、板、殼等）（整型） ，單元特性信息（厚度、內力矩等）（實型），材料信息（彈性模量，泊松比等）（實型）等。在FORTRAN 程序中，以上這些變量混合在一起，很難辨認，使程序的可讀性不好，如需要進行單元網絡的自適應劃分，節點

5、及單元的修改將非常困難。在進行C語言編譯過程中，采用結構struct 使每個節點信息存儲在一個結構體數組中，提高程序的可讀性，使數據結構更趨于合理。1.1. 八節點矩形單元介紹八節點矩形單元編號如Error! Reference source not found.所示圖 21八節點矩形單元的位移函數為： 其形函數為 式和式中，并且采用等參變換，則單元的坐標變換式可取為 單元應變矩陣為 式一般簡寫為 其中的子塊矩陣為 由于是、的函數，在中的、要按照復合函數來求導，即 從而有 因此，單元應力矩陣為 單元剛度矩陣為 其中積分采用三點高斯積分， （高斯積分點的總數），和或是加權系數，和是單元內的坐標.

6、。對于三點高斯積分，高斯積分點的位置： ，。單元等效節點荷載 結構剛度矩陣 結構結點荷載列陣 注意，對于式和式中的理解不是簡單的疊加而是集成?？倓偲胶夥匠?從式求出 將回代入式和式，得到和。1.2. 有限元分析的模塊組織一個典型的有限元分析過程主要包括以下幾個步驟：1) 讀輸入數據，定義節點及單元數組。2) 由邊界條件計算方程個數，賦值荷載列陣。3) 讀入在帶狀存儲的總剛度矩陣中單元和載荷信息。4) 定義總剛度陣數組。5) 組裝總剛度陣。6) 解方程組。輸入邊界條件(對稱條件)形成各荷載工況的節點荷載陣 總剛分解 回代求出位移及輸出 計算應變、應力形成單元剛陣單剛向總剛投放坐標變換輸入原始參數

7、計算總剛規模形成總剛方程向總節點荷載陣投放形成單元荷載陣調整幾何、彈性矩陣調整單元位移列陣圖 22其流程圖可見Error! Reference source not found.。3. 程序變量及函數說明1.3. 主要變量說明：主要程序變量說明 wide分析模型的寬 hight分析模型的高wsize寬度方向網格劃分尺寸hsize高度方向網格劃分尺寸npoin節點總數nelem單元總數struct node500節點結構體變量struct element 500單元結構體變量ifpre500節點約束信息posgp3高斯積分點weigp3權系數gpcod29高斯積分點總體坐標bmatx316單元變

8、形矩陣dmatx33單元彈性矩陣xjacm22雅可比矩陣xjaci22雅可比矩陣的逆矩陣djacb雅可比矩陣行列式的值estif136單元剛度矩陣shape8單元形函數deriv28形函數對局部坐標的導數cartd28形函數對整體坐標的導數A30000總體剛度矩陣eload16等效節點荷載gpcod29高斯積分點的總體坐標1.4. 主要函數說明主要函數說明void meshing（ ）網格劃分void coordinate( )節點整體坐標計算void input( )信息輸入void stif( )形成單元剛度矩陣void sfr2（）計算形函數對當前積分點及形函數對局部坐標的到數值void

9、 jacobi2( )形成雅可比矩陣void modps( )形成彈性矩陣void bmatps( )形成應變矩陣void load( )計算等效節點荷載void asem( )形成整體剛度矩陣和整體荷載列陣void solve( )求解整體方程void stre( )計算單元應力4. 程序流程圖任何一個有限元程序處理過程，都可以劃分為三個階段：1) 前處理：讀入數據和生成數據。2) 處理器：有限元的矩陣計算、組集和求解。3) 后處理：輸出節點位移、單元應變等。為了更清楚地描述程序，我們給出了程序的流程圖。5. 程序應用與ANSYS分析的比較為了驗證程序的正確性，我們取了一個簡單框架結構，分別

10、用自編程序和ANSYS進行分析和比較。1.5. 問題說明Error! Reference source not found.所示一簡支梁，高3m，長18m，承受均布荷載，取m，作為平面應力問題。圖 51圖 52由于結構和荷載對稱，只對右邊一半進行有限單元計算，如Error! Reference source not found.所示，而在y軸各節點布置水平連桿支座。1.6. ANSYS分析結果ANSYS分析中，采用plane82單元，在板單元上邊施加均布荷載，在y軸上的各結點約束UX方向的自由度，在板單元右下角的結點約束UX自由度，結點布置、網格劃分方案如Error! Reference so

11、urce not found.所示。圖 53Error! Reference source not found.和Error! Reference source not found.分別給出了結構的位移圖和應力云圖。圖 54 位移圖圖 55 各單元圖從Error! Reference source not found.和Error! Reference source not found.可以看到，跨中的位移和正應力最大，與簡支梁承受均布荷載，跨中撓度和正應力最大的力學常識相符合，可以初步認為模型是正確的。Error! Reference source not found.給出了的截面上的正應力

12、和Error! Reference source not found.給出了處各橫截面方向位移，其中 的單位為，的單位為。表格 1考查點的y/m1.51.00.50-0.5-1.0-1.5正應力-270.20-178.25-88.570.0088.57178.25270.21表格 2考查點的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3485-0.3380-0.3069-0.2571-0.1914-0.113301.7. 自編程序分析結果用自編程序進行分析時，采用了整個結構模型進行計算，其他條件不變，因編程水平有限，在后處理階段沒有給出節點處的位移與應力，而只能得到高

13、斯積分點處的位移與應力，得到積分點處的應力后，根據數值分析相關知識采用了插值進行處理，從而得到的截面上的正應力和 處各橫截面方向位移，分別列出表格如下：表格 3考查點的y/m1.51.00.50-0.5-1.0-1.5正應力-297.2-196.06-93.250.0093196.08297.23表格 4考查點的x/m01.53.04.56.07.59.0方向位移(10-6)-0.3481-0.3376-0.3065-0.2568-0.1910-0.112501.8. 結果分析比較為了更直觀的比較自編程序和ANSYS的計算結果，將Error! Reference source not foun

14、d.和Error! Reference source not found.的數據繪入Error! Reference source not found.，將Error! Reference source not found.和Error! Reference source not found.的數據繪入Error! Reference source not found.。圖 56 應力比較圖 57 位移比較 自編程序所得結果與ANSYS分析結果進行比較發現，應力偏大而位移偏小。分析其中的原因，一方面是編程水平有限，自編程序有很多不完善之處，有很多因素沒有考慮（溫度、變形、非線性等），所得結果可

15、信度不是很高，只能得到應力以及位移大概的分布情況；另一方面是在后處理階段，在對高斯積分點結果進行處理時，誤差難以避免，還有一方面的原因是在進行求解時保留數據精度不夠，造成誤差較大，并且進行數值運算時可能會造成大數吃小數，從而引起結果的差異。 參考文獻1 (美)史密斯(Smith,I.m.)著;王崧等譯.有限單元法編程(第三版)M.北京:電子工業出版社,20032 王勖成.有限單元法M.北京:清華大學出版社,20033 宋建輝,涂志剛.MATLAB語言及其在有限元編程中的應用J.湛江師范學院學報.2003.6(24):101-1054 鄭大玉, 連宏玉, 周躍發. 有限元編程與C語言J. 黑龍江

16、商學院學報.1997.3(13):23-285 王偉,劉德富.有限元編程中應用面向對象編程技術的探討J.三峽大學學報.2001.2(23):124-1286 汪文立, 呂士俊.二級C語言程序設計教程M. 北京:中國水利水電出版社,20067 趙翔龍.FORTRAN 90學習教程M.北京:北京大學出版社,2002附錄 程序源代碼#include<stdio.h>#include<math.h>/*The gemotry of the model*/ float mesh2; float wide,hight; float wsize,hsize,young,poiss,t

17、hick; int i,j,k,knelem,npoin,elem500,ielem; float coord21,props31,lnods8500,ifpre1; float posgp3,weigp3,estif136,elcod28; int npoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb; float gpcod29,bmatx316,dmatx33; float shape8,deriv28; float xjacm22,xjaci22,elcod28; float cartd28; float bmatx316,dmatx33; f

18、loat nload1; float press32,pgash2,dgash2;struct node int nodenu;/*The number of the node*/ float cor2;/*The coordinate of the node*/ float displacement2;/*The displacement of the node*/ float load2; /*The load of the node*/ float boundary2; node500; struct int elementnu;/*The number of element*/ int

19、 localnu8;/*Local number*/ int localcorx8; int localcory8; /*Local coordinate of node*/ float qx,qy,t;/*The stress and strain*/ element500;void meshing(float wide,float hight,float mesh2) printf("Please input the meshing sizen"); scanf("%f%f",&wsize,&hsize); getchar(); me

20、sh0=wide/wsize; mesh1=hight/hsize; /*The global coordinate of node*/ void coordinate() int i,j=1; float x,y; for(i=1;i<=2*mesh1+1;i+) x=0; while(x<=wide) if(i%2!=0) nodej.cor1+=wsize/2; nodej.cor2+=(i-1)*hsize; j+; else nodej.cor1+=wsize; nodej.cor2+=(i-1)*hsize; j+; main() int top500,bottom50

21、0;/*The number of top and bottom element*/ void input(); void stif(); void jacobi2( ); void load(); void asem(); void solve( ); void stre(); npoin=8; for(i=1;i<=8;i+) lnodsi1=i;for(i=1;i<=3;i+) scanf("%f",&propsi1); printf("The EX is %enThe uo is %8fnThe bt is %8f",prop

22、s11,props21,props31); getch(); printf("Please input the gemotry of the modeln"); scanf("%f%f",&wide,&hight); getchar(); printf("The wide is %0.3f,the hight is %0.3f",wide,hight); getchar(); meshing(wide,hight,mesh); printf("The wide size is %f,the hight siz

23、e is %fn",mesh0,mesh1); input(); getchar(); nelem=mesh0*mesh1; for(i=1;i<=nelem;i+) /*The element number*/ elementi.elementnu=i; for(i=1;i<mesh0;i+) npoin+=5; for(i=1;i<mesh1;i+) npoin+=3*mesh0+2; printf("%d",npoin); for(i=1;i<=npoin;i+) nodei.nodenu=i; for(i=1;i<=nelem

24、;i+) for(j=1;j<=8;j+) elementi.localnuj=j; for (i=1;i<=mesh0;i+) bottomi=i; j=1; for(i=mesh0*(mesh1-1)+1;i<=nelem;i+) topj+=i; for(i=1;i<j;i+) printf("the top number is %dn",topi); printf("%dn",element1.localnu8); coordinate(); stif(); jacobi2( ); load(); asem(); solv

25、e( ); stre();getchar();/*data input*/void input() int nnode=8; int notreadchar;/*input element node number*/printf("input element node numbern");for(ielem=1;ielem<=nelem;ielem+) for(i=1;i<=nnode;i+) scanf("%d",&lnodsiielem);/*Input restriction message*/printf("doub

26、le-digit is symbol of the restriction,1 representates stable,2 representates gived displacementn"); i=1; while(i) scanf("%d",&i); if(i!=0) scanf("%d",&j); scanf("%d",&ifprej); else break; /*Element stiffness matrix*/void stif() int kevab,nstre,ievab,ist

27、re,lnode,jstre,jevab; float kgasp,exisp,etasp, dvolu; float btdbm,dbmat3; void sfr();/*Gauss point*/ posgp1=-sqrt(0.6); posgp2=0; posgp3=sqrt(0.6); /*weight coefficient*/ weigp1=5.0/9.0; weigp2=8.0/9.0; weigp3=5.0/9.0; /*number of stress*/ nstre=3; /*form elastic matrix*/ for(ielem=1;ielem<=nelem

28、;ielem+) printf("The number of element is %dn",ielem);for(i=1;i<=8;i+) lnode=lnodsiielem; for(j=1;j<=2;j+) coordjlnode=nodelnode.corj; elcodji=coordjlnode; young=props11; poiss=props21; thick=props31; modps(young,poiss); /*The initial value of k*/ kevab=0; for(i=1;i<=16;i+) for(j=

29、1;j<=i;j+) kevab+; estifkevab=0.0; /*The gauss point shape function and deriv*/ kgasp=0; for(igaus=1;igaus<=3;igaus+) for(jgaus=1;jgaus<=3;jgaus+) kgasp=kgasp+1; exisp=posgpigaus; etasp=posgpjgaus; printf("The number of gauss point is %dn",kgasp); sfr2(exisp,etasp); jacob2(ielem,d

30、jacb,kgasp,elcod); dvolu=djacb*weigpigaus*weigpjgaus*thick; bmatps()/*The initial value of elastic matrix D*/ kevab=0; for(ievab=1;ievab<=16;ievab+) for(istre=1;istre<=nstre;istre+) dbmatistre=0.0; for(jstre=1;jstre<=nstre;jstre+) dbmatistre=dbmatistre+bmatxjstreievab*dmatxjstreistre; for(j

31、evab=1;jevab<=ievab;ievab+) kevab=kevab+1; btdbm=0.0; for(istre=1;istre<=nstre;istre+) btdbm=btdbm+dbmatistre*bmatxistrejevab; estifkevab=estifkevab+btdbm*dvolu; /*float sfr2(float s,float t) float s2,t2,ss,tt,st,stt,sst,st2;/*Shape function*/*these variables are to simplify the formula */ s2=

32、s*2.0; t2=t*2.0; ss=s*s; tt=t*t; st=s*t; stt=s*t*t; sst=s*s*t; st2=s*t*2.0; /*shape function*/ shape1=(-1.0+st+ss+tt-sst-stt)/4.0; shape2=(1.0-t-ss+sst)/2.0; shape3=(-1.0-st+ss+tt-sst+stt)/4.0; shape4=(1.0+s-tt-stt)/2.0; shape5=(-1.0+st+ss+tt+sst+stt)/4.0; shape6=(1.0+t-ss-sst)/2.0; shape7=(-1.0-st+

33、ss+tt+sst-stt)/4.0; shape8=(1.0-s-tt+stt)/2.0; /* differential coefficient to local coordinate*/ deriv11=(t+s2-st2-tt)/4.0; deriv12=-s+st; deriv13=(-t+s2-st2+tt)/4.0; deriv14=(1.0-tt)/2.0; deriv15=(t+s2+st2+tt)/4.0; deriv16=-s-st; deriv17=(-t+s2+st2-tt)/4.0; deriv18=(-1.0+tt)/2.0; deriv21=(s+t2-ss-s

34、t2)/4.0; deriv22=(-1.0+ss)/2.0; deriv23=(-s+t2-ss+st2)/4.0; deriv24=-t-st; deriv25=(s+t2+ss+st2)/4.0; deriv26=(1.0-ss)/2.0; deriv27=(-s+t2+ss-st2)/4.0; deriv28=-t+st; */*Jacobi matrix*/void jacobi2( ) /* xjacm22:jacobi matrix;xjaci22:j-1;elcod28:local coordinate*/ float cartd28,gpcod29; int i,j,k; f

35、or(i=1;i<=2;i+) gpcodikgasp=0.0; for(j=1;j<=8;j+) gpcodikgasp=gpcodikgasp+elcodij*shapej; for(i=1;i<=2;i+) for(j=1;j<=2;j+) xjacmij=0.0; for(k=1;k<=8;k+) xjacmij=xjacmij+derivik*elcodjk; /*Caculate the value of Jacobi matrix*/ djacb=xjacm11*xjacm22-xjacm12*xjacm21; printf("The de

36、t of jacobi matrix is %fn",djacb); if(djacb<=1e-6) printf("The jacobi det of %d is less or equal 0n",ielem); /*The element of j-1*/ xjaci11=xjacm22/djacb; xjaci22=xjacm11/djacb; xjaci12=-xjacm12/djacb; xjaci21=-xjacm21/djacb; /*The deriv of shape function to global coordinate*/ for

37、(i=1;i<=2;i+) for(k=1;k<=8;k+) cartdik=0.0; for(j=1;j<=2;j+) cartdik=cartdik+xjaciij*derivjk; /* Form elastic matris */*float modps() float bmatx316,dmatx33; float constant; int i,j,y,p; y=props11; p=props21; for(i=1;i<=3;i+) for(j=1;j<=3;j+) dmatxij=0.0; /*If Ntype=1,it is plane stre

38、ss question*/constant=y/(1.0-p*p); dmatx11=constant; dmatx22=constant; dmatx12=constant*p; dmatx21=constant*p; dmatx33=constant*(1.0-p)/2.0; */*Form strain matrix*/*void bmatps() float cartd28,shape8,deriv28,bmatx316,dmatx33; int npoin,nelem,ngash,mgash; float gpcod29;/*The initial value ofB*/ for(i

39、=1;i<=16;i+) for(j=1;j<=3;j+) bmatxji=0.0; /*Form B*/ ngash=0; for(i=1;i<=8;i+) mgash=ngash+1; ngash=mgash+1; bmatx1mgash=cartd1i; bmatx1ngash=0.0; bmatx2mgash=0.0; bmatx2ngash=cartd2i; bmatx3mgash=cartd2i; bmatx3ngash=cartd1i; */ /*equal load of node*/void load()float nload500; float elcod

40、23,eload16,posgp3,weigp3; float press32,pgash2,dgash2,noprs3; int iedge,nedge,kount; float sgmar,tau,s,dvolu,pxcom,pycom,n,m,t; int lnode,number,nloca; /*The number of load side*/ /*element load*/ /*The position of gauss point*/ printf("input the number of load side*/n"); scanf("%d&qu

41、ot;,&number); posgp1=-sqrt(0.6); posgp2=0.0; posgp3=sqrt(0.6); /*The weight of gauss point*/ weigp1=5.0/9.0; weigp2=8.0/9.0; weigp3=5.0/9.0; /*The initial load of element*/ for(i=1;i<=nelem;i+) nloadi=0.0; for(i=1;i<=number;i+) /*The initial value of equal node load*/ for(i=1;i<=16;i+)

42、eloadi=0; scanf("%d",&nedge); for(iedge=1;iedge<=nedge;iedge+) for(i=1;i<=3;i+) scanf("%f%f",&sgmar,&tau); /*if sgmar=0,input node load q(1,2,3)and t(1,2,3)*/ if(fabs(sgmar)<1e-8)&&(fabs(tau)<1e-8) for(j=1;j<=3;j+) for(i=1;i<=2;i+) scanf(&quo

43、t;%f",&pressji); else for(i=1;i<=3;i+) pressi1=sgmar; pressi2=tau; /*The coordinate of the load node*/ for(i=1;i<=3;i+) /*input load node*/ printf("input node load numbern"); scanf("%d",&noprsi); lnode=noprsi; for(j=1;j<=2;j+) coordjlnode=nodelnode.corj; elc

44、odji=coordjlnode; t=-1.0; for(igaus=1;igaus<=3;igaus+) s=posgpigaus; sfr(s,t); for(j=1;j<=2;j+) pgashj=0.0; dgashj=0.0; /*current point q,t and the value of deriv*/ pgashj=pgashj+pressij*shapei; dgashj=dgashj+elcodji*deriv1i; dvolu=weigpigaus; pxcom=dgash1*pgash2-dgash2*pgash1; pycom=dgash1*pg

45、ash1+dgash2*pgash2; for(i=1;i<=8;i+) nloca=lnodsiielem; if(nloca=noprs1) j=i+2; break; j=i+2; kount=0; for(k=i;k<=j;k+) kount=kount+1; n=(k-1)*2+1; m=n+1; /*if the local node numner >8,then it is 1*/ if(k>8) n=1; m=2; /*sum to get the equal node load*/ eloadn=eloadn+shapekount*pxcom*dvolu; eloadm=eloadm+shapekount*pycom*dvolu; prin