Fortran平面钢架有限元分析报告.docx
- 文档编号:5895492
- 上传时间:2023-01-02
- 格式:DOCX
- 页数:21
- 大小:229.20KB
Fortran平面钢架有限元分析报告.docx
《Fortran平面钢架有限元分析报告.docx》由会员分享,可在线阅读,更多相关《Fortran平面钢架有限元分析报告.docx(21页珍藏版)》请在冰豆网上搜索。
Fortran平面钢架有限元分析报告
1有限元分析软件的开发
1.1程序功能
该程序为平面刚架静力分析程序,能针对平面刚架间问题进行有限元计算,计算杆端位移及杆端力大小。
程序从磁盘文件中读取单元编号、节点编号及坐标、材料属性、荷载、边界条件等信息;将杆端位移,杆端力等计算结果以磁盘文件的形式输出,采用等带宽二维数组存储整体刚度矩阵并使用高斯消去法进行求解。
1.2程序结构及流程
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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Fortran 平面 钢架 有限元分析 报告