Fortran平面钢架有限元分析.docx
《Fortran平面钢架有限元分析.docx》由会员分享，可在线阅读，更多相关《Fortran平面钢架有限元分析.docx（22页珍藏版）》请在冰豆网上搜索。
Fortran平面钢架有限元分析
1有限元分析软件的开发
1.1程序功能
该程序为平面刚架静力分析程序，能针对平面刚架间问题进行有限元计算，计算杆端位移及杆端力大小。
程序从磁盘文件中读取单元编号、节点编号及坐标、材料属性、荷载、边界条件等信息；将杆端位移，杆端力等计算结果以磁盘文件的形式输出，采用等带宽二维数组存储整体刚度矩阵并使用高斯消去法进行求解。
1.2程序结构及流程
开始
标题及数组说明
（读入题目序号NO）
NO是否为零?
形成整体刚度矩阵
计算并打印各杆轴力
解方程并打印杆端位移
形成结点载荷
结束
子程序READ
子程序MKE
子程序MAKE
子程序MR
子程序MF
是
否
读入数据并打印
子程序MULV6
子程序CALM
子程序MK
子程序MULV
子程序TARN
子程序SOLV
子程序MADE
子程序PE
1.3程序的输入与输出
详细介绍输入输出数据的格式。
如：
数据文件分几个部分，各有几行，分别包含哪些容及其类型、先后次序，等等。
输入，共有九行。
第一行：
7,13,5,1,2,2。
分别为，7个结点，13个自由度，5个单元，1个类型，2个结点荷载，2个非结点荷载。
第二行：
1,2,3,0.0,0.0,0,0,,6.0,0.0。
分别为：
一号结点的位移序号，x方向为1，y方向为2，转角为3，坐标为（0.0,0.0），因为二号结点固结在地面，所以二号结点的位移序号，x方向为0，y方向为0，转角为0，坐标为（6.0,0.0）。
第三行：
4,5,6,0.0,6.0,4,5,7,0.0,6.0。
分别为：
三号结点的位移序号，x方向为4，y方向为5，转角为6，坐标为（0.0,6.0）,四号结点位移序号x方向和y相同，转角为7，坐标为（,0.0,6.0）。
第四行：
8,9,10,6.0,6.0,0,0,11,0.0,12.0.五号结点位移序号，x方向为8，y方向为9，转角为10，坐标为（6.0,6.0）。
因为六号结点铰接在地面，所以六号结点的位移序号，x方向和y方向为0，转角为11，坐标为（0.0,12.0）。
第五行：
12,0,13,6.0,12.0.因为七号结点与地面用滑动支座固定，所以七号结点的位移序号，x方向为12，y方向为0，转角为13，坐标（6.0,12.0）.
第六行：
1,2,1,1,3,1,4,5,1,3,6,1,5,7,1,分别为，1号和2号结点组成的单元为1号类型。
1号和3号结点组成的单元为1号类型，4号和5号结点组成的为1号类型，3号和6号结点组成的单元为1号类型，5号和7号结点组成的单元为1号类型。
第七行：
分别为，弹性模量为E＝2×108kN/m2，截面面积A=0.16m2,惯性矩I=0.002m4。
第八行：
1号结点转角方向的集中力偶为-20.0kN，3号结点集中力为10.0KN。
第九行：
1号单元，受集中力（集中力型号为3），大小为15.0kN，到始端的距离为3.0。
5号单元，受均布力（均布力型号为1），大小为5.0kN，到端点的距离为5.0。
第十行：
0为计算终止符。
输出：
第一部分为输入的数据。
RESULTSOFCALCULATION以下为输出结果，第二部分的第一段为4个结点的x,y方向的位移和转角。
第二段为1,2,3号单元的轴力，剪力和弯矩。
1.4程序求解中遇到的问题
1对实例进行计算时，坐标原点选用不同的点，会导致整个题目的坐标值发生改变，输入的容会有所不同，最后的结果也不相同
2对结点荷载和非结点荷载的正负判断不同，结点荷载的方向和整体坐标有关，非结点荷载方向判断和局部坐标有关。
3在非结点荷载中，均布荷载和集中力到始端的距离判断不同。
2有限元分析算例
2.1算例说明
已知图示刚架，各杆的材料及截面均相同，弹性模量E＝2×108kN/m2,A=0.16m2,惯性矩I=0.002m4，q=5kN/m.，一号单元集中力为15KN,一号结点集中力偶为20KN*M,三号结点集中力为10KN.试求刚架的力。
节点编号如图
2.2理论分析
对所选取的力学问题进行理论分析，要有详细的推导过程和计算结果。
1力计算
对结构进行分析，可以看出1,2,4单元组成的是二次超静定结构，3,5单元是静定结构。
因此先对3,5单元组成的结构进行分析。
如上图所示，可以根据x，y方向力平衡，对结点七力矩平衡算得支座反力。
再画出其弯矩，剪力轴力图。
然后对1,2,4单元组成的结构分析。
用力法解超静定，将结点六的约束解除，加上支座反力x1=1，x2=1.画出M1，M2,MP图。
MP图M1图M2图
然后画出其弯矩，剪力，轴力图
弯矩图剪力图轴力图
2位移计算
计算结点1位移，x方向加单位力1
其剪力与弯矩图都为零，轴力图为
根据公式:
轴力图图乘Δ1x=（12.92*6*1）/EA=（12.92*6*1）/3*10^7*0.16=1.615e-5
y方向加单位力1，忽略剪力的影响，弯矩图图乘
Δ1y=-1/EI（3*15.05*3/2+1/2*3*2.82*2/3*3+1/2*1.27*（3+1.27/3）*17.87）+1/EI（1/2*1.73*24.32*（4.27+2*1.73/3））=1.611*E-5（略小于程序结果）
加单位力偶，剪力与轴力图为零，弯矩图为
θ1=1/EI（3*15.05*1+1/2*3*2.82*1+1/2*1.27*1*17.87）-1/EI（1/2*1.73*24.32*1）=0.6617e-3
计算结点3的位移，x方向加单位力1.
弯矩轴力
Δ1X=-1/EI（3*15.05*6+1/2*3*2.82*6+1/2*1.27*6*17.87）+1/EI（1/2*1.73*24.32*6）-1/EI（1/2*2.17*35.05*（3.29+2*2.71/3））+1/EI（1/2*3.29*42.47*2.71/3）+（12.92*6*1）/EA=-0.675E-2.其它位移同理可得。
2.3输入输出数据
输入：
输出：
2.4分析结果
理论分析中，因为力计算应用了力法，所以程序所得结果和理论结果一致。
而对位移进行理论分析时忽略了剪力的影响，所以理论位移略小于程序所计算的结果。
可以看出软件的的正确性很高，但是此软件只适用计算平面杆系结构，不能解决弹性力学问题.
结点1位移
u
v
θ
理论值
1.615e-5
1.611*E-5
0.6617e-3
程序值
1.615e-5
1.640*E-5
0.6617e-3
3程序源代码
附上完整的程序源代码。
PROGRAMPFAP
CANALYSISPROGRAMFORPLANEFRAME
REALK（200,200）,KE（6,6）,AKE（6,6）,X（100）,Y（100）,AL（100）,
#EAI（3,100）,PJ（100）,PF（2,100）,R（6,6）,P（100）,FF（6）,
#FE（6）,D（100）,ADE（6）,DE（6）,RT（6,6）,AFE（6）,F（3）
INTEGERJE（2,100）,JN（3,100）,JPJ（100）,JPF（2,100）,M（6）,JEAI（100）,NO
OPEN（6,FILE='ht2.TXT'）
OPEN（8,FILE='ht.txt',STATUS='NEW'）
1READ（6,*）NO
IF（NO.EQ.0）STOP
WRITE（8,'（/9X,A5,I3,A1）'）'（NO=',NO,'）'
CALLREAD（NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,JPJ,PJ,JPF,PF）
DO5I=1,N
P（I）=0.0
DO5J=1,N
5K（I,J）=0.0
DO10IE=1,NEL
CALLMKE（KE,IE,JE,JEAI,EAI,X,Y,AL）
CALLMR（R,IE,JE,X,Y）
CALLMAKE（KE,R,AKE）
CALLCALM（M,IE,JN,JE）
CALLMK（K,AKE,M）
10CONTINUE
DO20IP=1,NPF
CALLMR（R,JPF（1,IP）,JE,X,Y）
CALLTRAN（R,RT）
CALLPE（FE,IP,JPF,PF,AL）
CALLMULV6（RT,FE,AFE）
CALLCALM（M,JPF（1,IP）,JN,JE）
CALLMF（P,AFE,M）
20CONTINUE
DO30I=1,NPJ
30P（JPJ（I））=P（JPJ（I））+PJ（I）
CALLSOLV（K,P,D,N）
WRITE（8,'（/2（26（1H*）,A））'）'RESULTSOFCALCULATION'
WRITE（8,'（/28X,A）'）'NODELDISPLACEMENT'
WRITE（8,40）
40FORMAT（9X,'NO.N',4X,'X-DISPLACEMENT',2X,
#'Y-DISPLACEMENT',3X,'ANG.ROT.（RAD）'）
DO60KK=1,NJ
DO50II=1,3
F（II）=0.0
I1=JN（II,KK）
50IF（I1.GT.0）F（II）=D（I1）
60WRITE（8,70）KK,F
（1）,F
（2）,F（3）
70FORMAT（4X,I8,2X,3G16.5）
WRITE（8,'（/30X,A）'）'ELEMANTFORCES'
WRITE（8,80）
80FORMAT（2X,'NO.E',4X,'N
（1）',9X,'Q
（1）',9X,'M
（1）',
#9X,'N
（2）',9X,'Q
（2）',9X,'M
（2）'）
DO120IE=1,NEL
CALLMADE（IE,JN,JE,D,ADE）
CALLMKE（KE,IE,JE,JEAI,EAI,X,Y,AL）
CALLMR（R,IE,JE,X,Y）
CALLMULV6（R,ADE,DE）
CALLMULV6（KE,DE,FF）
DO100IP=1,NPF
IF（JPF（1,IP）.EQ.IE）THEN
CALLPE（FE,IP,JPF,PF,AL）
DO90I=1,6
90FF（I）=FF（I）-FE（I）
ENDIF
100CONTINUE
WRITE（8,110）IE,（FF（I）,I=1,6）
110FORMAT（I5,2X,6G13.7）
120CONTINUE
GOTO1
END
SUBROUTINEREAD（NJ,N,NEL,NM,NPJ,NPF,JN,X,Y,JE,JEAI,EAI,
#JPJ,PJ,JPF,PF）
REALX（100）,Y（100）,EAI（3,100）,PJ（100）,PF（2,100）
INTEGERJE（2,100）,JN（3,100）,JPJ（100）,JPF（2,100）,JEAI（100）,
#TITLE（20）
READ（6,'（20A4）'）（TITLE（I）,I=1,20）
WRITE（8,'（/9X,20A4）'）TITLE
READ（6,*）NJ,N,NEL,NM,NPJ,NPF
WRITE（8,'（/3（8X,A4,1H:
I2））'）'NJ=',NJ,
#'N=',N,'NE=',NEL,'NM=',NM,'NPJ=',NPJ,'NPF=',NPF
WRITE（8,10）
10FORMAT（/8X,'NO.N
（1）
（2）（3）',10X,'X',9X,'Y'）
READ（6,*）（（JN（J,I）,J=1,3）,X（I）,Y（I）,I=1,NJ）
DO20I=1,NJ
20WRITE（8,'（8X,1H（,I2,1H）,3I6,4X,2F10.3）'）I,JN（1,I）,JN（2,I）,
#JN（3,I）,X（I）,Y（I）
READ（6,*）（JE（1,I）,JE（2,I）,JEAI（I）,I=1,NEL）
WRITE（8,25）
25FORMAT（/7X,'NO.E
（1）
（2）NO.MATNO.E
（1）
（2）NO.MAT'）
N2=（NEL+1）/2
DO30I=1,N2-1
30WRITE（8,40）I,（JE（J,I）,J=1,2）,JEAI（I）,
#I+N2,（JE（J,I+N2）,J=1,2）,JEAI（I+N2）
IF（N2*2.NE.NEL）WRITE（8,40）N2,JE（1,N2）,JE（2,N2）,JEAI（N2）
IF（N2*2.EQ.NEL）WRITE（8,40）N2,JE（1,N2）,JE（2,N2）,JEAI（N2）,
#NEL,JE（1,NEL）,JE（2,NEL）,JEAI（NEL）
40FORMAT（4X,4I6,2X,4I6）
READ（6,*）（（EAI（I,J）,I=1,3）,J=1,NM）
WRITE（8,50）（J,（EAI（I,J）,I=1,3）,J=1,NM）
50FORMAT（/3X,'NO.MAT',6X,'ELASTICMODULUS',5X,
#'AREA',7X,'MOMENTOFINERTIA'/（I6,5X,3G16.4））
IF（NPJ.EQ.0）GOTO70
WRITE（8,'（/20X,12HNODELLOADS）'）
WRITE（8,'（16XA）'）'NO.DISP.VALUE'
READ（6,*）（JPJ（I）,PJ（I）,I=1,NPJ）
DO60I=1,NPJ
60WRITE（8,'（14X,I7,F16.3）'）JPJ（I）,PJ（I）
70CONTINUE
IF（NPF.EQ.0）GOTO100
WRITE（8,'（/20X,16HNON-NODELLOADS）'）
WRITE（8,'（7X,A,8X,A,9X,A）'）'NO.ENO.LOAD.MODEL','A','C'
READ（6,*）（JPF（1,I）,JPF（2,I）,PF（1,I）,PF（2,I）,I=1,NPF）
DO80I=1,NPF
80WRITE（8,90）（JPF（J,I）,J=1,2）,PF（1,I）,PF（2,I）
90FORMAT（6X,I3,8X,I4,8X,2F10.3）
100CONTINUE
RETURN
END
SUBROUTINEMKE（KE,IE,JE,JEAI,EAI,X,Y,AL）
REALKE（6,6）,X（100）,Y（100）,EAI（3,100）,AL（100）,L
INTEGERJE（2,100）,JEAI（100）
II=JE（1,IE）
JJ=JE（2,IE）
MT=JEAI（IE）
L=SQRT（（X（JJ）-X（II））**2+（Y（JJ）-Y（II））**2）
AL（IE）=L
A1=EAI（1,MT）*EAI（2,MT）/L
A2=EAI（1,MT）*EAI（3,MT）/L**3
A3=EAI（1,MT）*EAI（3,MT）/L**2
A4=EAI（1,MT）*EAI（3,MT）/L
KE（1,1）=A1
KE（1,4）=-A1
KE（2,2）=12*A2
KE（2,3）=6*A3
KE（2,5）=-12*A2
KE（2,6）=6*A3
KE（3,3）=4*A4
KE（3,5）=-6*A3
KE（3,6）=2*A4
KE（4,4）=A1
KE（5,5）=12*A2
KE（5,6）=-6*A3
KE（6,6）=4*A4
DO10I=1,6
DO10K=1,6
10KE（K,I）=KE（I,K）
RETURN
END
SUBROUTINEMR（R,IE,JE,X,Y）
REALR（6,6）,X（100）,Y（100）,L,CX,CY
INTEGERJE（2,100）
I=JE（1,IE）
J=JE（2,IE）
L=SQRT（（X（J）-X（I））**2+（Y（J）-Y（I））**2）
CX=（X（J）-X（I））/L
CY=（Y（J）-Y（I））/L
DO10J=1,6
DO10I=1,6
10R（I,J）=0.0
DO20I=1,4,3
R（I,I）=CX
R（I,I+1）=CY
R（I+1,I）=-CY
R（I+1,I+1）=CX
20R（I+2,I+2）=1
RETURN
END
SUBROUTINEMAKE（KE,R,AKE）
REALKE（6,6）,R（6,6）,RT（6,6）,TMP（6,6）,AKE（6,6）
CALLTRAN（R,RT）
CALLMULV（RT,KE,TMP）
CALLMULV（TMP,R,AKE）
RETURN
END
SUBROUTINECALM（M,IE,JN,JE）
INTEGERM（6）,JN（3,100）,JE（2,100）,IE
DO10I=1,3
M（I）=JN（I,JE（1,IE））
10M（I+3）=JN（I,JE（2,IE））
RETURN
END
SUBROUTINEMK（K,AKE,M）
REALK（200,200）,AKE（6,6）
INTEGERM（6）
DO10I=1,6
DO10J=1,6
IF（M（I）.NE.0.AND.M（J）.NE.0）
#K（M（I）,M（J））=K（M（I）,M（J））+AKE（I,J）
10CONTINUE
RETURN
END
SUBROUTINEPE（FE,IP,JPF,PF,AL）
REALFE（6）,PF（2,100）,AL（100）,L
INTEGERJPF（2,100）
A=PF（1,IP）
C=PF（2,IP）
L=AL（JPF（1,IP））
IND=JPF（2,IP）
DO5I=1,6
5FE（I）=0.0
GOTO（10,20,30,40,50,60）,IND
10FE
（2）=（7.*A/20.+3.*C/20.）*L
FE（3）=（A/20.+C/30.）*L**2
FE（5）=（3.*A/20.+7*C/20.）*L
FE（6）=-（A/30.+C/20.）*L**2
RETURN
20FE（5）=A*C**3*（2.*L-C）/2./L**3
FE
（2）=A*C-FE（5）
FE（3）=A*C**2*（6.*L*L-8.*C*L+3.*C*C）/12./L/L
FE（6）=-A*C**3*（4.*L-3.*C）/12./L/L
RETURN
30FE
（2）=A*（L-C）**2*（L+2.*C）/L**3
FE（3）=A*C*（L-C）**2/L**2
FE（5）=A-FE
（2）
FE（6）=-A*C**2*（L-C）/L**2
RETURN
40FE
（2）=-6.*A*C*（L-C）/L**3
FE（3）=A*（L-C）*（L-3.*C）/L**2
FE（5）=-FE
（2）
FE（6）=A*C*（3.*C-2.*L）/L**2
RETURN
50FE
（1）=A*（1.-C/L）
FE（4）=A*C/L
RETURN
60FE
（1）=C*L/2.
FE（4）=FE
（1）
RETURN
END
SUBROUTINEMULV6（A,B,C）
REALC（6）,A（6,6）,B（6）
DO10I=1,6
C（I）=0.0
DO10J=1,6
10C（I）=C（I）+A（I,J）*B（J）
RETURN
END
SUBROUTINEMF（P,AFE,M）
REALP（100）,AFE（6）
INTEGERM（6）
DO10I=1,6
IF（M（I）.NE.0）P（M（I））=AFE（I）+P（M（I））
10CONTINUE
RETURN
END
SUBROUTINESOLV（AK,P,D,N）
REALAK（200,200）,P（100）,D（100）
DO5I=1,100
5D（I）=P（I）
DO10K=1,N-1
DO10I=K+1,N
C=-AK（K,I）/AK（K,K）
DO20J=I,N
20AK（I,J）=AK（I,J）+C*AK（K,J）
10D（I）=D（I）+C*D（K）
D（N）=D（N）/AK（N,N）
DO40I=N-1,1,-1
DO30J=I+1,N
30D（I）=D（I）-AK（I,J）*D（J）
40D（I）=D（I）/AK（I,I）
RETURN
END
SUBROUTINEMADE（IE,JN,JE,D,ADE）
REALADE（6）,D（100）
INTEGERIE,JN（3,100）,JE（2,100）
DO3I=1,6
3ADE（I）=0.0
DO10I=1,3
IF（JN（I,JE（1,IE））.NE.0）ADE（I）=D（JN（I,JE（1,IE）））
IF（JN（I,JE（2,IE））.NE.0）ADE（I+3）=D（JN（I,JE（2,IE）））
10CONTINUE
RETURN
END
SUBROUTINETRAN（R,RT）
REALR（6,6）,RT（6,6）
DO10I=1,6
DO10J=1,6
10RT（I,J）=R（J,I）
RETURN
END
SUBROUTINEMULV（A,B,C）
REALA（6,6）,B（6,6）,C（6,6）
DO10I=1,6
DO10J=1,6
C（I,J）=0.0
DO10K=1,6
10C（I,J）=C（I,J）+A（I,K）*B（K,J）
RETURN
END
升级会员 