有限元程序设计.doc

.平面问题有限元程序设计理学院 学号 兔兔设计人 完成日期 一、 程序功能说明本程序适用于节点荷载作用下的桁架分析问题，当有节间荷载存在时可按照静力等效原理将其转化为节点荷载。可求解平面桁架在静力荷载作用下的内力和位移。二、 框图的设计开 始 输 入 数 据 数 组 定 义计算各杆截面面积和半带宽 调用形成单刚矩阵unit调用形成半带宽存贮的结构原始刚度矩阵total有节点荷载否输入节点荷载值，并将其送入相应的荷载列阵p（n）中考虑结构是否自重将杆自重引起的等效荷载叠加到p（n）中支座处理、解方程，并输出u（n）、v（n）调用unit，求各单元杆端内力 结 束单元循环没有有否是 三、 程序的标识符及数组说明npoin 最大节点数nelem 最大单元数 nload 节点的荷载总数nzero 节点的约束位移总数wt 结构的自重ee 材料的弹性模量ll 一维数组，用于存放单元杆件的长度aa 一维数组，用于存放单元杆件的面积coord 节点坐标数组lnode 单元节点数组bh 二维数组，用于存放单元截面尺寸nres 二维数组，用于存放约束的位移值jp 二维数组，用于存放节点的荷载值estif 四维数组，用于存放整体坐标系下的单元刚度矩阵astif 二十维数组，用于存放半带宽结构原始刚度矩阵p 用于存放节点的荷载列阵u 用于存放节点x方向的位移值v 用于存放节点y方向的位移值四、 源程序integer 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 输入已知数据data npoin,nelem,njp,nres,ee,wt/3,3,1,3,21000,0/data coord/0,6,0,0,0,6/data lnode/1,2,1,2,3,3/data bh/3*2,3*10/data res/3*0,1,2,4/c 计算各单元面积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计算半带宽 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=2*(nhbw+1) write(2,220) nhbw220 format(1x,nhbw=, i2)c单元循环 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=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,荷载总数,8x,水平荷载,8x,铅垂荷载) do 730 k=1,no730 write(2,740) 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,nhbw-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)900 continue write(2,910)910 format(/10x,节点位移,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,单元号,8x,节点号,8x,n(kn),8x,q(kn)/) do 980 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,i2,8x,f8.4,8x,f8.4/18x,i2,8x,f8.4,8x,f8.4)980 continue write(2,1000)1000format(/28x,结束,/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-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(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)=-estif(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，不计各杆的自重，求在荷载作用下，各杆的轴力。图1桁架例图计算模型的输入数据有以下各量：npoinnresnelemnloadyogwt3331210000.0节点单元12lnode(1,i)12lnode(2,i)23lnode(3,i)13编号数组123res(i,1)0.00.00.0res(i,2)124节点 i123coord(i,1)0.06.00.0coord(i，2)0.00.06.0jp(i,1)10.0jp(i,2)5单元号i123bh(i,1)101010bh(i,2)222输出