有限元程序設(shè)計(jì)

1、平面問題有限元程序設(shè)計(jì)理學(xué)院 學(xué)號(hào) 設(shè)計(jì)人 完成日期 一、 程序功能說明本程序適用于節(jié)點(diǎn)荷載作用下的桁架分析問題，當(dāng)有節(jié)間荷載存在時(shí)可按照靜力等效原理將其轉(zhuǎn)化為節(jié)點(diǎn)荷載?？汕蠼馄矫骅旒茉陟o力荷載作用下的內(nèi)力和位移。二、 框圖的設(shè)計(jì)開 始 輸 入 數(shù) 據(jù) 數(shù) 組 定 義計(jì)算各桿截面面積和半帶寬 調(diào)用形成單剛矩陣unit調(diào)用形成半帶寬存貯的結(jié)構(gòu)原始剛度矩陣total有節(jié)點(diǎn)荷載否輸入節(jié)點(diǎn)荷載值，并將其送入相應(yīng)的荷載列陣p（n）中考慮結(jié)構(gòu)是否自重將桿自重引起的等效荷載疊加到p（n）中支座處理、解方程，并輸出u（n）、v（n）調(diào)用unit，求各單元桿端內(nèi)力 結(jié) 束單元循環(huán)沒有有否是 三、 程序的標(biāo)識(shí)符及

2、數(shù)組說明npoin 最大節(jié)點(diǎn)數(shù)nelem 最大單元數(shù) nload 節(jié)點(diǎn)的荷載總數(shù)nzero 節(jié)點(diǎn)的約束位移總數(shù)wt 結(jié)構(gòu)的自重ee 材料的彈性模量ll 一維數(shù)組，用于存放單元桿件的長度aa 一維數(shù)組，用于存放單元桿件的面積coord 節(jié)點(diǎn)坐標(biāo)數(shù)組lnode 單元節(jié)點(diǎn)數(shù)組bh 二維數(shù)組，用于存放單元截面尺寸nres 二維數(shù)組，用于存放約束的位移值jp 二維數(shù)組，用于存放節(jié)點(diǎn)的荷載值estif 四維數(shù)組，用于存放整體坐標(biāo)系下的單元?jiǎng)偠染仃嘺stif 二十維數(shù)組，用于存放半帶寬結(jié)構(gòu)原始剛度矩陣p 用于存放節(jié)點(diǎn)的荷載列陣u 用于存放節(jié)點(diǎn)x方向的位移值v 用于存放節(jié)點(diǎn)y方向的位移值四、 源程序integ

3、er e,nelem,z,h real ll,estif,astif,jp dimension coord(3,2),lnode(3,2),aa(200),bh(3,2),res(3,2), &ll(200),estif(4,4),astif(400,20),jp(1,2),p(400),u(200), &v(200) open(2,file='d:nmxjia.dat',status='new') c 輸入已知數(shù)據(jù)data npoin,nelem,njp,nres,ee,wt/3,3,1,3,21000,0/data coord/0,6,0,0

4、,0,6/data lnode/1,2,1,2,3,3/data bh/3*2,3*10/data res/3*0,1,2,4/c 計(jì)算各單元面積do 200 e=1,nelemaa(e)=bh(e,1)*bh(e,2) call unit(e,ee,coord,lnode,aa,estif,ll,cx,cy)200 continue c計(jì)算半帶寬 l2=2*npoin nhbw=0 do 210 e=1,nelem m=abs(lnode(e,1)-lnode(e,2) if(nhbw.lt.m) nhbw=m210 continuewrite(2,*) '半帶寬' nhbw

5、=2*(nhbw+1) write(2,220) nhbw220 format(1x,'nhbw=', i2)c單元循環(huán) do 300 i1=1,l2 do 300 j1=1,nhbw300 astif(i1,j1)=0.0 do 400 e=1,nelem call unit(e,ee,coord,lnode,aa,estif,ll,cx,cy) call total(e,lnode,estif,astif) 400 continue do 560 n=1,l2560 p(n)=0.0 if(njp.eq.0) goto 650data jp/10,5/ do 630 k1=

6、1,njp nn=jp(k1,2)+0.1 630 p(nn)=jp(k1,1) 650 if(wt.le.0.0) goto 750 do 700 e=1,nelem n1=lnode(e,1) n2=lnode(e,2) p(2*n1)=p(2*n1)-wt*aa(e)*ll(e)/2.0 p(2*n2)=p(2*n2)-wt*aa(e)*ll(e)/2.0700 continue write(2,710)710 format(/4x,'荷載總數(shù)',8x,'水平荷載',8x,'鉛垂荷載') do 730 k=1,no730 write(2,7

7、40) k,p(2*k-1),p(2*k)740 format(4x,i2,8x,f8.3,8x,f8.3) 750 do 800 i1=1,nres z=res(i1,2)+1e-5 astif(z,1)=astif(z,1)*1e8 p(z)=astif(z,1)*res(i1,1)800 continue do 850 k1=1,l2-1 if(l2.gt.(k1+nhbw-1) then im=k1+nhbw-1 else im=l2 endif do 850 i1=k1+1,im l1=i1-k1+1 c1=astif(k1,l1)/astif(k1,1) do 830 j1=1,n

8、hbw-l1+1 mm=j1+i1-k1 astif(i1,j1)=astif(i1,j1)-c1*astif(k1,mm)830 continue p(i1)=p(i1)-c1*p(k1)850 continue p(l2)=p(l2)/astif(l2,1) do 900 i1=l2-1,1,-1 if(nhbw.gt.(l2-i1+1) then jm=l2-i1+1 else jm=nhbw endif do 880 j1=2,jm h=j1+i1-1 p(i1)=p(i1)-astif(i1,j1)*p(h)880 continue p(i1)=p(i1)/astif(i1,1)90

9、0 continue write(2,910)910 format(/10x,'節(jié)點(diǎn)位移',10x,'水平位移',10x,'鉛垂位移'/) do 930 n=1,no u(n)=p(2*n-1) v(n)=p(2*n)930 write(2,950) n,u(n),v(n)950 format(15x,i2,6x,f12.7,6x,f12.7) write(2,970)970 format(/4x,'單元號(hào)',8x,'節(jié)點(diǎn)號(hào)',8x,'n(kn)',8x,'q(kn)'/) do 9

10、80 e=1,nelem call unit(e,ee,coord,lnode,aa,estif,ll,cx,cy) n1=lnode(e,1) n2=lnode(e,2) ulnode=u(n1)-u(n2) vlnode=v(n1)-v(n2) d1=estif(1,1)*ulnode+estif(1,2)*vlnode d2=estif(1,2)*ulnode+estif(2,2)*vlnode fi=cx*d1+cy*d2 fj=-fi ti=-cy*d1+cx*d2 tj=-ti write(2,990) e,n1,fi,ti,n2,fj,tj990 format(4x,i2,12x

11、,i2,8x,f8.4,8x,f8.4/18x,i2,8x,f8.4,8x,f8.4)980 continue write(2,1000)1000format(/28x,'結(jié)束',/15x,35('*')/) stop end 子程序total形成總剛度矩陣subroutine total(e,lnode,estif,astif) integer e,dh,zl,dl real estif,astif dimension lnode(3,2),estif(4,4),astif(400,20) do 40 i1=1,2 do 40 ii=1,2 kh=2*(i1-

12、1)+ii dh=2*(lnode(e,i1)-1)+ii do 40 j1=1,2 do 40 jj=1,2 kl=2*(j1-1)+jj zl=2*(lnode(e,j1)-1)+jj dl=zl-dh+1 if(dl.gt.0) astif(dh,dl)=astif(dh,dl)+estif(kh,kl)40 continue return end 子程序unit形成單剛subroutine unit(e,ee,coord,lnode,aa,estif,ll,cx,cy) integer e real ll,estif dimension coord(3,2),lnode(3,2),aa

13、(200),ll(200),estif(4,4) n1=lnode(e,1) n2=lnode(e,2) cx=coord(n2,1)-coord(n1,1) cy=coord(n2,2)-coord(n1,2) ll(e)=sqrt(cx*cx+cy*cy) cx=cx/ll(e) cy=cy/ll(e) eal=ee*aa(e)/ll(e) estif(1,1)=eal*cx*cx estif(1,2)=eal*cx*cy estif(2,2)=eal*cy*cy estif(2,1)=estif(1,2) do 10 i=1,2 do 10 j=1,2 estif(i,j+2)=-est

14、if(i,j) estif(i+2,j)=-estif(i,j)10 estif(i+2,j+2)=estif(i,j) return end五、 算例如圖1所示桁架，已知桿件材料的彈性模量,桿件截面高度h=10cm，截面的寬度為b=2cm，不計(jì)各桿的自重，求在荷載作用下，各桿的軸力。圖1桁架例圖計(jì)算模型的輸入數(shù)據(jù)有以下各量：npoinnresnelemnloadyogwt3331210000.0節(jié)點(diǎn)單元12lnode(1,i)12lnode(2,i)23lnode(3,i)13編號(hào)數(shù)組123res(i,1)0.00.00.0res(i,2)124節(jié)點(diǎn) i123coord(i,1)0.06.00.0coord(i，2)0.00.06.0jp(i,1)10.0jp(i,2)5單元號(hào)i123bh(i,1)101010bh(