Fortran语言有限元程序分析报告平面钢架.docx
《Fortran语言有限元程序分析报告平面钢架.docx》由会员分享，可在线阅读，更多相关《Fortran语言有限元程序分析报告平面钢架.docx（16页珍藏版）》请在冰豆网上搜索。
Fortran语言有限元程序分析报告平面钢架
程序框图：
程序特点：
问题类型：
可用于计算结构力学的平面刚架问题
单元类型：
直接利用杆单元
载荷类型：
节点载荷及非节点载荷，其中非节点载荷包括均布荷载和垂直于杆件的集中荷载
材料性质：
所有杆单元几何性质相同，且由相同的均匀材料组成
方程求解：
结构刚度矩阵采用满阵存放，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）存放荷载方向代码（1—x方向；2—y方向；3—转角），pj（np,3）存放荷载大小
pf（ne,4）:
存放非结点荷载信息，其中pf（ne,1）存放荷载作用单元号，pf（ne,2）存放荷载代码（1—均布荷载，2—垂直集中荷载），pf（ne,3）存放荷载大小，pf（ne,4）荷载作用距离（均布荷载，集中荷载均以单元起始结点为计算起始位置）。
2.子例行子程序哑元信息：
第一部分：
基本部分
I.subroutinelsc（Length&Sin&Cos）:
输入哑元：
m（单元号）,nj,ne,x,y,ij
输出哑元：
bl（杆件长度）,si（正弦值）,co（余弦值）
II.subroutineelv（ElementLocationVector）:
输入哑元：
m,ne,nj,ij,jn
输出哑元：
lv（单元定位数组）
III.subroutineesm（ElementStiffnessMatrix）:
输入哑元：
e,a,zi,bl,si,co
输出哑元：
ek（整体坐标系下的单刚矩阵）
IV.subroutineeff（ElementFixed-endForces）
输入哑元：
i,pf,nf,bl
输出哑元：
fo（局部坐标系下单元固端力）
第二部分：
主程序直接调用部分
I.subroutinetsm（TotalStiffnessMatrix计算总刚矩阵）
输入哑元：
ne,nj,n,e,x,y,ij,a,zi,jn
输出哑元：
tk
II.subroutinejlp（JointLoadVector计算结点荷载）
输入哑元：
ne,nj,n,np,nf,x,y,ij,jn,pj,pf
输出哑元：
p（结点荷载列矩阵）
III.subroutinegauss（带列主元素消去的高斯法）
输入（输出）哑元：
tk,p,n;（注意，算出位移后，直接存储到结点荷载列矩阵）
IV.subroutinemvn（Member-endforcesofelements计算各单元的杆端力）
输入哑元：
ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p
3.文件管理：
源程序文件：
pff.for
程序需读入的数据文件：
input.txt
程序输出的数据文件：
output
4.数据文件格式：
【输入文件格式】：
栏目
格式说明
实际需输入的数据
基本模型数据
第1行，每两个数之间用“，”号隔开
单元个数，结点个数，总自由度，弹性模量，杆截面积，惯性矩，结点荷载个数，非结点荷载个数
结点位置信息
第2行，每两个数之间用“，”号隔开
依次输入各结点的坐标（x,y）
单元结点信息
每输入一个单元换行（回车），两个数之间用“，”号隔开
依次输入各单元的起点结点号和终点结点号
结点约束信息
每输入一个结点换行（回车），两个数之间用“，”号隔开
按先处理法要求，输入各结点编号
结点荷载信息
每个结点荷载成一行，每两个数之间用“，”号隔开
每行依次输入荷载作用的结点号，荷载方向代码，荷载大小（参考“主要变量”的叙述）
非结点荷载信息
每个非结点荷载成一行，每两个数之间用“，”号隔开
每行依次输入荷载作用单元，荷载代码，荷载大小，荷载作用”长度”（参考“主要向量”的叙述）
【输出文件格式】：
1.第1部分：
每行数据依次为：
结点号，结点x方向位移，结点y方向位移，结点转角位移
2.第2部分：
每行数据依次为：
单元号，
，
，
，
，
，
源程序：
programPFF
implicitnone
realtk（100,100）,x（50）,y（50）,p（100）,pj（50,3）,pf（50,4）
integerij（50,2）,jn（50,3）
integerne,nj,n,np,nf
reale,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
callinput（ne,nj,x,y,ij,jn,np,nf,pj,pf）
calltsm（ne,nj,n,e,x,y,ij,a,zi,jn,tk）
calljlp（ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p）
callgauss（tk,p,n）
callmvn（ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p）
end
subroutineinput（ne,nj,x,y,ij,jn,np,nf,pj,pf）
dimensionx（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（np>0）read（1,*）（（pj（i,j）,j=1,3）,i=1,np）
if（nf>0）read（1,*）（（pf（i,j）,j=1,4）,i=1,nf）
end
subroutinetsm（ne,nj,n,e,x,y,ij,a,zi,jn,tk）
dimensionx（nj）,y（nj）,ij（ne,2）,jn（nj,3）,tk（n,n）,ek（6,6）,lv（6）
doi=1,n
doj=1,n
tk（i,j）=0
enddo
enddo
dom=1,ne
calllsc（m,ne,nj,x,y,ij,bl,si,co）
callesm（e,a,zi,bl,si,co,ek）
callelv（m,ne,nj,ij,jn,lv）
dol=1,6
i=lv（l）
if（i/=0）then
dok=1,6
j=lv（k）
if（j/=0）tk（i,j）=tk（i,j）+ek（l,k）
enddo
endif
enddo
enddo
end
subroutinelsc（m,ne,nj,x,y,ij,bl,si,co）
dimensionx（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/bl
end
subroutineesm（e,a,zi,bl,si,co,ek）
dimensionek（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=c2
ek（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*s6
doi=1,5
doj=i+1,6
ek（j,i）=ek（i,j）
enddo
enddo
end
subroutineelv（m,ne,nj,ij,jn,lv）
dimensionij（ne,2）,jn（nj,3）,lv（6）
i=ij（m,1）
j=ij（m,2）
dok=1,3
lv（k）=jn（i,k）
lv（k+3）=jn（j,k）
enddo
end
subroutinejlp（ne,nj,n,np,nf,x,y,ij,jn,pj,pf,p）
dimensionx（nj）,y（nj）,ij（ne,2）,jn（nj,3）,pj（np,3）,pf（nf,4）,p（n）,fo（6）,pe（6）,lv（6）
doi=1,n
p（i）=0.0
enddo
if（np>0）then
doi=1,np
j=int（pj（i,1））
k=int（pj（i,2））
l=jn（j,k）
if（l/=0）p（l）=pj（i,3）
enddo
endif
if（nf>0）then
doi=1,nf
m=int（pf（i,1））
calllsc（m,ne,nj,x,y,ij,bl,si,co）
calleff（i,pf,nf,bl,fo）
callelv（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）
doj=1,6
l=lv（j）
if（l/=0）p（l）=p（l）+pe（j）
enddo
enddo
endif
end
subroutineeff（i,pf,nf,bl,fo）
dimensionpf（nf,4）,fo（6）
no=int（pf（i,2））
q=pf（i,3）
c=pf（i,4）
b=bl-c
c1=c/bl
c2=c1*c1
c3=c1*c2
doj=1,6
fo（j）=0.0
enddo
goto（10,20）,no
10fo
（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）
return
20fo
（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*b
return
end
subroutinegauss（e,d,n）
dimensione（n,n）,d（n）,a（n,n+1）
doi=1,n
doj=1,n
a（i,j）=e（i,j）
enddo
enddo
doi=1,n
a（i,n+1）=d（i）
enddo
dok=1,n-1
doi=k+1,n
if（abs（a（i,k））>abs（a（k,k）））then
doj=1,n+1
c=a（k,j）
a（k,j）=a（i,j）
a（i,j）=c
enddo
else
endif
enddo
doi=k+1,n
a（i,k）=a（i,k）/a（k,k）
doj=k+1,n+1
a（i,j）=a（i,j）-a（i,k）*a（k,j）
enddo
enddo
enddo
a（n,n+1）=a（n,n+1）/a（n,n）
doi=n-1,1,-1
doj=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=0
enddo
doi=1,n
d（i）=a（i,n+1）
enddo
end
subroutinemvn（ne,nj,n,nf,e,x,y,ij,a,zi,jn,pf,p）
dimensionx（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）
10format（//2x,"结点位移"/5x,"结点号",9x,"u向位移",9x,"v向位移",9x,"角位移"）
doj=1,nj
doi=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）
20format（2x,i6,4x,3e15.6）
enddo
write（2,30）
30format（/2x,"单元杆端力及弯矩"/4x,"单元号",13x,"Fx",17x,"Fy",17x,"弯矩"）
dom=1,ne
calllsc（m,ne,nj,x,y,ij,bl,si,co）
callesm（e,a,zi,bl,si,co,ek）
callelv（m,ne,nj,ij,jn,lv）
doi=1,6
l=lv（i）
d（i）=0.0
if（l/=0）d（i）=p（l）
enddo
doi=1,6
fd（i）=0.0
doj=1,6
fd（i）=fd（i）+ek（i,j）*d（j）
enddo
enddo
f
（1）=fd
（1）*co+fd
（2）*si
f
（2）=-fd
（1）*si+fd
（2）*co
f（3）=fd（3）
f（4）=fd（4）*co+fd（5）*si
f（5）=-fd（4）*si+fd（5）*co
f（6）=fd（6）
if（nf>0）then
doi=1,nf
l=int（pf（i,1））
if（m==l）then
calleff（i,pf,nf,bl,fo）
doj=1,6
f（j）=f（j）+fo（j）
enddo
endif
enddo
endif
write（2,40）m,f
40format（2x,i8,4x,"Ix=",f12.4,3x,"Iy=",f12.4,3x,"Mi=",f12.4/14x,"Jx=",f12.4,3x,"Jy=",f12.4,3x,"Mj=",f12.4）
enddo
end
【算例】：
课题二：
平面刚架有限元程序分析
题目一：
分析如图所示结构，其中
，
，
，
，
。
用程序PFF计算此结构的输入数据文件如下：
输出的计算结果文件如下：
【学习体会】：
1.通过这次课程设计，我进一步了解了关于Fortran90的一些更复杂更深入的用法，例如子例行子程序的调用，主程序与子程序之间数据的传递等。
2.了解了有限元设计的模块化（结构化）设计。
结构化程序设计事程序结构清晰，设计工作可分步进行。
3.进一步了解了杆系结构的有限元设计。
知道了程序处理方法和人工处理方法的不同之处。
【参考书目】：
《土木工程结构分析程序设计》赵更新中国水利水电出版社2002
《有限单元法》丁科北京大学出版社2006
升级会员 