Fortran语言 有限元程序分析 平面钢架.doc

程序框图：Sub-lsc 单元长及角度Sub-elv 单元定位向量Sub-esm 单刚矩阵Sub-lsc 单元长及角度Sub-elv 单元定位向量Sub-eff 单元固端力Sub-lsc 单元长及角度Sub-elv 单元定位向量Sub-esm 单刚矩阵Sub-eff 单元固端力开始调用input输入数据tsm 形成结构刚度矩阵jlp 形成结点荷载列阵gauss 解方程组mvn 计算和输出杆端力结束输出output数据文件程序特点：问题类型：可用于计算结构力学的平面刚架问题单元类型：直接利用杆单元载荷类型：节点载荷及非节点载荷，其中非节点载荷包括均布荷载和垂直于杆件的集中荷载材料性质：所有杆单元几何性质相同，且由相同的均匀材料组成方程求解：结构刚度矩阵采用满阵存放，Gauss消元过程采用数值分析中的列主元素消去法输入文件：按先处理法的要求，由手工生成输入数据文件1.主要变量：ne: 单元个数nj: 结点个数n: 自由度e: 弹性模量（单位：KN/m2）a: 杆截面积zi: 惯性矩np: 结点荷载个数nf: 非结点荷载个数x(nj): 存放结点的x轴坐标y(nj): 存放结点的y轴坐标ij(ne,2): 存放单元结点编号，其中ij(nj,1)存放起始结点编号，ij(nj,2)存放终止结点编号jn(nj,3): 存放结点位移编号，以组成单元定位数组pj(np,3): 存放结点荷载信息，其中pj(np,1)存放结点荷载作用结点号，pj(np,2)存放荷载方向代码(1x方向；2y方向；3转角)，pj(np,3)存放荷载大小pf(ne,4): 存放非结点荷载信息，其中pf(ne,1)存放荷载作用单元号，pf(ne,2)存放荷载代码(1均布荷载，2垂直集中荷载)，pf(ne,3)存放荷载大小，pf(ne,4)荷载作用距离（均布荷载，集中荷载均以单元起始结点为计算起始位置）。2.子例行子程序哑元信息： 第一部分：基本部分I. subroutine lsc(Length & Sin & Cos): 输入哑元：m(单元号),nj,ne,x,y,ij 输出哑元：bl(杆件长度),si(正弦值),co(余弦值) II. subroutine elv(Element Location Vector): 输入哑元：m,ne,nj,ij,jn输出哑元：lv(单元定位数组)III. subroutine esm(Element Stiffness Matrix): 输入哑元：e,a,zi,bl,si,co输出哑元：ek(整体坐标系下的单刚矩阵)IV. subroutine eff(Element Fixed-end Forces) 输入哑元：i,pf,nf,bl输出哑元：fo(局部坐标系下单元固端力) 第二部分：主程序直接调用部分I. subroutine tsm(Total Stiffness Matrix 计算总刚矩阵)输入哑元：ne,nj,n,e,x,y,ij,a,zi,jn输出哑元：tkII. subroutine jlp(Joint Load Vector 计算结点荷载) 输入哑元：ne,nj,n,np,nf,x,y,ij,jn,pj,pf输出哑元：p(结点荷载列矩阵)III. subroutine gauss(带列主元素消去的高斯法)输入(输出)哑元：tk,p,n ;(注意，算出位移后，直接存储到结点荷载列矩阵)IV. subroutine mvn(Member-end forces of elements 计算各单元的杆端力) 输入哑元：ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p3.文件管理： 源程序文件：pff.for 程序需读入的数据文件：input.txt 程序输出的数据文件：output4.数据文件格式： 【输入文件格式】：栏目格式说明实际需输入的数据基本模型数据第1行，每两个数之间用“，”号隔开单元个数，结点个数，总自由度，弹性模量，杆截面积，惯性矩，结点荷载个数，非结点荷载个数结点位置信息第2行，每两个数之间用“，”号隔开依次输入各结点的坐标（x,y）单元结点信息每输入一个单元换行（回车），两个数之间用“，”号隔开依次输入各单元的起点结点号和终点结点号结点约束信息每输入一个结点换行（回车），两个数之间用“，”号隔开按先处理法要求，输入各结点编号结点荷载信息每个结点荷载成一行，每两个数之间用“，”号隔开每行依次输入荷载作用的结点号，荷载方向代码，荷载大小（参考“主要变量”的叙述）非结点荷载信息每个非结点荷载成一行，每两个数之间用“，”号隔开每行依次输入荷载作用单元，荷载代码，荷载大小，荷载作用”长度”（参考“主要向量”的叙述）【输出文件格式】：1. 第1部分： 每行数据依次为：结点号，结点x方向位移，结点y方向位移，结点转角位移2. 第2部分： 每行数据依次为：单元号，源程序：program PFF implicit none real tk(100,100),x(50),y(50),p(100),pj(50,3),pf(50,4) integer ij(50,2),jn(50,3) integer ne,nj,n,np,nf real e,a,zi open(1,file=input.txt,status=old) open(2,file=output.txt,status=old) read(1,*) ne,nj,n,e,a,zi,np,nf call input(ne,nj,x,y,ij,jn,np,nf,pj,pf) call tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk) call jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p) call gauss(tk,p,n) call mvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p)endsubroutine input(ne,nj,x,y,ij,jn,np,nf,pj,pf)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pj(np,3),pf(nf,4)read(1,*)(x(i),y(i),i=1,nj)read(1,*)(ij(i,1),ij(i,2),i=1,ne)read(1,*)(jn(i,j),j=1,3),i=1,nj)if (np0) read(1,*)(pj(i,j),j=1,3),i=1,np)if (nf0) read(1,*)(pf(i,j),j=1,4),i=1,nf)endsubroutine tsm(ne,nj,n,e,x,y,ij,a,zi,jn,tk)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),tk(n,n),ek(6,6),lv(6)do i=1,n do j=1,n tk(i,j)=0 enddoenddodo m=1,ne call lsc(m,ne,nj,x,y,ij,bl,si,co) call esm(e,a,zi,bl,si,co,ek) call elv(m,ne,nj,ij,jn,lv) do l=1,6 i=lv(l) if (i/=0) then do k=1,6 j=lv(k) if (j/=0) tk(i,j)=tk(i,j)+ek(l,k) enddo endifenddoenddoendsubroutine lsc(m,ne,nj,x,y,ij,bl,si,co) dimension x(nj),y(nj),ij(ne,2) i=ij(m,1) j=ij(m,2) dx=x(j)-x(i) dy=y(j)-y(i) bl=sqrt(dx*dx+dy*dy) si=dy/bl co=dx/blendsubroutine esm(e,a,zi,bl,si,co,ek) dimension ek(6,6) c1=e*a/bl c2=2.0*e*zi/bl c3=3.0*c2/bl c4=2.0*c3/bl s1=c1*co*co+c4*si*si s2=(c1-c4)*si*co s3=c3*si s4=c1*si*si+c4*co*co s5=c3*co s6=c2ek(1,1)=s1 ek(1,2)=s2 ek(1,3)=-s3 ek(1,4)=-s1 ek(1,5)=-s2 ek(1,6)=-s3 ek(2,2)=s4 ek(2,3)=s5 ek(2,4)=-s2 ek(2,5)=-s4 ek(2,6)=s5 ek(3,3)=2*s6 ek(3,4)=s3 ek(3,5)=-s5 ek(3,6)=s6 ek(4,4)=s1 ek(4,5)=s2 ek(4,6)=s3 ek(5,5)=s4 ek(5,6)=-s5 ek(6,6)=2.0*s6do i=1,5 do j=i+1,6 ek(j,i)=ek(i,j) enddo enddoendsubroutine elv(m,ne,nj,ij,jn,lv) dimension ij(ne,2),jn(nj,3),lv(6) i=ij(m,1) j=ij(m,2) do k=1,3 lv(k)=jn(i,k) lv(k+3)=jn(j,k) enddoendsubroutine jlp(ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pj(np,3),pf(nf,4),p(n),fo(6),pe(6),lv(6)do i=1,n p(i)=0.0enddoif (np0) then do i=1,np j=int(pj(i,1) k=int(pj(i,2) l=jn(j,k) if (l/=0) p(l)=pj(i,3) enddoendifif(nf0) then do i=1,nf m=int(pf(i,1) call lsc(m,ne,nj,x,y,ij,bl,si,co) call eff(i,pf,nf,bl,fo) call elv(m,ne,nj,ij,jn,lv) pe(1)=-fo(1)*co+fo(2)*si pe(2)=-fo(1)*si-fo(2)*co pe(3)=-fo(3) pe(4)=-fo(4)*co+fo(5)*si pe(5)=-fo(4)*si-fo(5)*co pe(6)=-fo(6) do j=1,6 l=lv(j) if (l/=0) p(l)=p(l)+pe(j) enddoenddo endifendsubroutine eff(i,pf,nf,bl,fo)dimension pf(nf,4),fo(6)no=int(pf(i,2)q=pf(i,3)c=pf(i,4)b=bl-cc1=c/blc2=c1*c1c3=c1*c2do j=1,6 fo(j)=0.0enddogoto(10,20),no10 fo(2)=-q*c*(1.0-c2+c3/2.0) fo(3)=-q*c*c*(0.5-2.0*c1/3.0+0.25*c2) fo(5)=-q*c*c2*(1.0-0.5*c1) fo(6)=q*c*c*c1*(1.0/3.0-0.25*c1)return20 fo(2)=-q*b*b*(1.0+2.0*c1)/bl/bl fo(3)=-q*c*b*b/bl/bl fo(5)=-q*c2*(1.0+2.0*b/bl) fo(6)=q*c2*breturnendsubroutine gauss(e,d,n)dimension e(n,n),d(n),a(n,n+1)do i=1,n do j=1,n a(i,j)=e(i,j) enddoenddodo i=1,n a(i,n+1)=d(i)enddodo k=1,n-1 do i=k+1,n if (abs(a(i,k)abs(a(k,k) then do j=1,n+1 c=a(k,j) a(k,j)=a(i,j) a(i,j)=c enddo else endif enddo do i=k+1,n a(i,k)=a(i,k)/a(k,k) do j=k+1,n+1 a(i,j)=a(i,j)-a(i,k)*a(k,j) enddo enddoenddoa(n,n+1)=a(n,n+1)/a(n,n)do i=n-1,1,-1 do j=i+1,n p=p+a(i,j)*a(j,n+1) enddo a(i,n+1)=(a(i,n+1)-p)/a(i,i) p=0enddodo i=1,n d(i)=a(i,n+1)enddoendsubroutine mvn(ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p)dimension x(nj),y(nj),ij(ne,2),jn(nj,3),pf(nf,4),lv(6),fo(6),d(6),fd(6),f(6),ek(6,6),p(n)write(2,10)10 format(/2x,结点位移/5x,结点号,9x,u向位移,9x,v向位移,9x,角位移)do j=1,nj do i=1,3 d(i)=0.0 l=jn(j,i) if (l/=0) d(i)=p(l) enddo write(2,20)j,d(1),d(2),d(3)20 format(2x,i6,4x,3e15.6)enddowrite(2,30)30 format(/2x,单元杆端力及弯矩/4x,单元号,13x,Fx,17x,Fy,17x,弯矩)do m=1,ne call lsc(m,ne,nj,x,y,ij,bl,si,co) call esm(e,a,zi,bl,si,co,ek) call elv(m,ne,nj,ij,jn,lv) do i=1,6 l=lv(i) d(i)=0.0 if(l/=0) d(i)=p(l) enddo do i=1,6 fd(i)=0.0 do j=1,6 fd(i)=fd(i)+ek(i,j)*d(j) enddo enddo f(1)