Fortran语言编写的有限元结构程序

1、 算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。各杆ea，ei相同。已知： 计算简图如下：（1）输入原始数据控制参数 3,5,8,7,1,2（ne,nj,n,nw,npj,npf）结点坐标集结点未知量编号 单元杆端结点编号及单元ea、ei 结点荷载 非结点荷载 （2）输出结果 ne= 3 nj= 5 n= 8 nw= 7 npj= 1 npf= 2node x y xx yy zz 1 0.0000 0.0000 0 0 0 2 0.0000 4.0000 1 2 3 3 0.0000 4.0000 1 2 4 4 4.0000 4.0000 5 6 7 5 4.0000 0.

2、0000 0 0 8 element node-i node-j ea ei 1 1 2 0.400000e+07 0.160000e+05 2 3 4 0.400000e+07 0.160000e+05 3 5 4 0.400000e+07 0.160000e+05 code px-py-pm 7. -15.0000 element ind a q 1. 2. 2.0000 -18.0000 2. 1. 4.0000 -25.0000 node u v ceta 1 0.000000e+00 0.000000e+00 0.000000e+00 2 -0.221743e-02 -0.46461

3、9e-04 -0.139404e-02 3 -0.221743e-02 -0.464619e-04 0.357876e-02 4 -0.222472e-02 -0.535381e-04 -0.298554e-02 5 0.000000e+00 0.000000e+00 0.658499e-03 element n q m 1 n1= 46.4619 q1= 10.7119 m1= -6.8477 n2= -46.4619 q2= 7.2881 m2= 0.0000 2 n1= 7.2881 q1= 46.4619 m1= 0.0000 n2= -7.2881 q2= 53.5381 m2= 1

4、4.1523 3 n1= 53.5381 q1= 7.2881 m1= 0.0000 n2= -53.5381 q2= -7.2881 m2= -29.1523算例二计算简图及结果输出用平面刚架静力计算程序下图结构的内力。已知：桁架单元的抗拉刚度为，平面刚架单元的抗拉刚度为已知：，抗弯刚度为。计算简图如下：（1）输入原始数据控制参数 （ne,nj,n,nw,npj,npf）结点坐标集结点未知量编号 单元杆端结点编号及单元ea、ei 非结点荷载 （2）输出结果 ne= 5 nj= 4 n= 8 nw= 7 npj= 0 npf= 1 node x y xx yy zz 1 0.0000 0.00

5、00 0 0 1 2 4.0000 0.0000 2 3 4 3 4.0000 -3.0000 5 6 0 4 8.0000 0.0000 7 0 8 element node-i node-j ea ei 1 1 2 0.600000e+07 0.184000e+06 2 2 4 0.600000e+07 0.184000e+06 3 3 1 0.200000e+07 0.000000e+00 4 3 2 0.200000e+07 0.000000e+00 5 3 4 0.200000e+07 0.000000e+00 element ind a q 1. 1. 4.0000 -20.000

6、0 node u v ceta 1 0.000000e+00 0.000000e+00 0.312593e-03 2 -0.202759e-04 -0.253871e-03 -0.144928e-03 3 -0.202759e-04 -0.185440e-03 0.000000e+00 4 -0.405518e-04 0.000000e+00 -0.227378e-04 element n q m 1 n1= 30.4138 q1= 37.1896 m1= 0.0000 n2= -30.4138 q2= 42.8104 m2= 11.2415 2 n1= 30.4138 q1= 2.8104

7、m1= -11.2415 n2= -30.4138 q2= -2.8104 m2= 0.0000 3 n1= -38.0173 q1= 0.0000 m1= 0.0000 n2= 38.0173 q2= 0.0000 m2= 0.0000 4 n1= 45.6207 q1= 0.0000 m1= 0.0000 n2= -45.6207 q2= 0.0000 m2= 0.0000 5 n1= -38.0173 q1= 0.0000 m1= 0.0000 n2= 38.0173 q2= 0.0000 m2= 0.0000c 主程序c (一)输入原始数据 dimension je(2,100),jn

8、(3,100),jc(6),ea(100),ei(100),x(100), $y(100),pj(2,50),pf(4,100) real*8 ke(6,6),kd(6,6),t(6,6),p(300),kb(200,20),f(6),fo(6), $d(6),bl,si,co,s,c open (5,file='rpf1.txt')open (6,file='jieguo1.dat',status='new')read (5,*) ne,nj,n,nw,npj,npfread (5,*) (x(j),y(j),(jn(i,j),i=1,3),j

9、=1,nj)read (5,*) (je(i,j),i=1,2),ea(j),ei(j),j=1,ne)if (npj.ne.0) read (5,*) (pj(i,j),i=1,2),j=1,npj)if (npf.ne.0) read (5,*) (pf(i,j),i=1,4),j=1,npf)write (6,10) ne,nj,n,nw,npj,npfwrite (6,20) (j,x(j),y(j),(jn(i,j),i=1,3),j=1,nj)write (6,30) (j,(je(i,j),i=1,2),ea(j),ei(j),j=1,ne)if (npj.ne.0) write

10、 (6,40) (pj(i,j),i=1,2),j=1,npj)if (npf.ne.0) write (6,50) (pf(i,j),i=1,4),j=1,npf)10 format(/6x,'ne=',i5,2x,'nj=',i5,2x,'n=',i5,2x,'nw=',i5,2x, $'npj=',i5,2x,'npf='i5)20 format(/7x,'node',7x,'x',11x,'y',12x,'xx',8x,'

11、;yy',8x,'zz'/ $(1x,i10,2f12.4,3i10)30format(/4x,'element',4x,'node-i',4x,'node-j',11x,'ea',13x,'ei'/ $(1x,3i10,2e15.6)40 format(/7x,'code',7x,'px-py-pm'/(1x,f10.0,f15.4)50format(/4x,'element',7x,'ind',10x,'a'

12、,14x,'q',/ $(1x,2f10.0,2f15.4)c (二)形成总结点荷载向量 do 55 i=1,n55 p(i)=0.00 if(npj.eq.0) go to 65do 60 i=1,npjl=pj(1,i)60 p(l)=pj(2,i)65 if(npf.eq.0) go to 90 do 70 i=1,npfm=pf(1,i)call scl(m,ne,nj,bl,si,co,je,x,y)call efx(i,npf,bl,pf,fo)call ctm(si,co,t)call ejc(m,ne,nj,je,jn,jc)do 75 l=1,6s=0.00d

13、o 80 k=1,680 s=s-t(k,l)*fo(k) f(l)=s75 continue do 85 j=1,6l=jc(j)if(l.eq.0) go to 85p(l)=p(l)+f(j)85 continue70 continuec (三)形成整体刚度矩阵90 do 95 i=1,n do 100 j=1,nw100 kb(i,j)=0.0095 continue do 105 m=1,necall scl(m,ne,nj,bl,si,co,je,x,y)call ctm(si,co,t)call esm(m,ne,bl,ea,ei,kd)call ejc(m,ne,nj,je,j

14、n,jc)do 110 i=1,6do 115 j=1,6s=0.00do 120 l=1,6do 125 k=1,6125 s=s+t(l,i)*kd(l,k)*t(k,j)120 continue ke(i,j)=s115 continue110 continue do 130 l=1,6i=jc(l)if(i.eq.0) go to 130do 135 k=1,6j=jc(k)if(j.eq.0.or.j.lt.i)go to 135jj=j-i+1kb(i,jj)=kb(i,jj)+ke(l,k)135 continue130 continue105 continuec (四)解线性方

15、程组 n1=n-1do 140 k=1,n1im=k+nw-1if(n.lt.im) im=ni1=k+1do 145 i=i1,iml=i-k+1c=kb(k,l)/kb(k,1)jm=nw-l+1do 150 j=1,jmjj=j+i-k150 kb(i,j)=kb(i,j)-c*kb(k,jj)145 p(i)=p(i)-c*p(k)140 continue p(n)=p(n)/kb(n,1)do 155 k=1,n1i=n-kjm=k+1if(nw.lt.jm) jm=nwdo 160 j=2,jml=j+i-1160 p(i)=p(i)-kb(i,j)*p(l)155 p(i)=p(

16、i)/kb(i,1) write(6,165)165 format(/7x,'node',10x,'u',14x,'v',11x,'ceta') do 170 i=1,njdo 175 j=1,3d(j)=0.00l=jn(j,i)if(l.eq.0) go to 175d(j)=p(l)175 continue write(6,180) i,d(1),d(2),d(3)180 format(1x,i10,3e15.6)170 continuec (五)求单元杆端内力 write(6,200)200 format(/4x,'

17、;element',13x,'n',17x,'q',17x,'m') do 205 m=1,necall scl(m,ne,nj,bl,si,co,je,x,y)call esm(m,ne,bl,ea,ei,kd)call ctm(si,co,t)call ejc(m,ne,nj,je,jn,jc)do 210 i=1,6l=jc(i)d(i)=0.00if(l.eq.0) go to 210d(i)=p(l)210 continue do 220 i=1,6f(i)=0.00do 230 j=1,6do 240 k=1,6240 f(i

18、)=f(i)+kd(i,j)*t(j,k)*d(k)230 continue220 continue if(npf.eq.0) go to 270do 250 i=1,npfl=pf(1,i)if(m.ne.l) go to 250call efx(i,npf,bl,pf,fo)do 260 j=1,6260 f(j)=f(j)+fo(j)250 continue270 write(6,280) m,(f(i),i=1,6)280 format(/1x,i10,3x,'n1=',f12.4,3x,'q1=',f12.4,3x,'m1=',f12.

19、4 $/14x,'n2=',f12.4,3x,'q2=',f12.4,3x,'m2=',f12.4)205continue close(5)stopendc 子程序c (六)形成单元定位向量 subroutine ejc(m,ne,nj,je,jn,jc)dimension je(2,ne),jn(3,nj),jc(6)j1=je(1,m)j2=je(2,m)do 10 i=1,3jc(i)=jn(i,j1)10 jc(i+3)=jn(i,j2) returnendc (七)求单元常数subroutine scl(m,ne,nj,bl,si,co

20、,je,x,y)dimension je(2,ne),x(nj),y(nj)real*8 bl,si,co,dx,dyj1=je(1,m)j2=je(2,m)dx=x(j2)-x(j1)dy=y(j2)-y(j1)bl=dsqrt(dx*dx+dy*dy)si=dy/blco=dx/blreturnendc (八)形成单元刚度矩阵 subroutine esm(m,ne,bl,ea,ei,kd)dimension ea(ne),ei(ne)real*8 kd(6,6),bl,s,g,g1,g2,g3g=ea(m)/blg1=2.00*ei(m)/blg2=3.00*g1/blg3=2.00*g

21、2/bldo 10 i=1,6do 10 j=1,610 kd(i,j)=0.00 kd(1,1)=gkd(1,4)=-gkd(4,4)=gkd(2,2)=g3kd(5,5)=g3kd(2,5)=-g3kd(2,3)=-g2kd(2,6)=-g2kd(3,5)=g2kd(5,6)=g2kd(3,3)=2.00*g1kd(6,6)=2.00*g1kd(3,6)=g1do 20 i=1,5 i1=i+1do 30 j=i1,630 kd(j,i)=kd(i,j)20 continue return endc (九)形成单元坐标转换矩阵 subroutine ctm(si,co,t)real*8 t(6,6),si,codo 10 i=1,6do 10 j=1,610 t(i,j)=0.00 t(1,1)=cot(1,2)=sit(2,1)=-sit(2,2)=cot(3,3)=1.00do 20 i=1,3do 20 j=1,320 t(i+3,j+3)=t(i,j) return endc (十)形成单元固端力 subroutine efx(i,npf,bl,pf,fo)dimension pf(4,npf)real*8 fo(6),a,b,c,g,q,s,blind=pf(2,i)a=pf(3,i)