有限元计算结构力学fortran程序.docx
- 文档编号:8737304
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:27
- 大小:26.15KB
有限元计算结构力学fortran程序.docx
《有限元计算结构力学fortran程序.docx》由会员分享,可在线阅读,更多相关《有限元计算结构力学fortran程序.docx(27页珍藏版)》请在冰豆网上搜索。
有限元计算结构力学fortran程序
有限元计算结构力学fortran程序
计算结构力学程序
计算结构力学编程大作业
时间,
2007年6月
!
!
!
****************************************************************************
!
!
!
关于程序的说明
!
!
!
****************************************************************************
!
一、功能:
!
1、可计算包括节点力,一般非节点力,支座沉降、温度荷载作用、制造误差的平
!
面桁架、梁、刚架及其组合结构的节点位移与杆端力;
!
2、可同时计算多种工况下的节点位移与杆端力。
!
*****************************************************************************
!
******************************************************************************
!
!
二、变量说明:
!
NE——单元数;
!
N——结构中自由度数;
!
NJ——节点数;
!
NS——特殊节点数,包括支座节点、主从节点(1节点不做主节点)、连接桁架的铰节点(没有转角);
!
NAI——结构的单元截面类型数;
!
MT——单元截面类型号;
!
NL——荷载工况数;
!
H——截面高度;
!
E——弹性模量;
!
JC——单元定位向量数组;
!
X(NJ),Y(NJ)——节点的X,Y坐标值;
!
JE(NE,2)——单元两端节点码数组;
!
AI(NAI,2)——按单元类型顺序存放A与I,AI(I,1)—第I类单元的截面积,AI(I,2)—第I类单元的
!
惯性矩;
!
MT(NE)——单元所属单元类型号;
!
JS(NS,4)——特殊节点信息,JS(I,1)—结点码;JS(I,2),JS(I,3),JS(I,4)—U,V,CETA约束信息,
!
有约束为1,没有约束为0;从节点某位移同主节点时位移时,该位移约束信息填主节点码;
!
!
PJ(NP,3)——节点荷载信息数组;PJ(I,1)—节点力所在节点号;PJ(I,2)—节点力作用坐标方向:
!
坐标方向U,V,M分别为1,2,3;PJ(I,3)—节点力的大小(含正负号);U,V方向集中力时,
!
与坐标轴正向同向为正,M按右手法则为正;本程序推导过程取y轴向下为正。
!
!
PF(NF,4)——非节点荷载数组,并给出以下类型说明:
!
前6类型数据输法(梯形等可以用叠加法计算):
!
PF(I,1)-单元码;PF(I,2)-类型;PF(I,3)-荷载大小;PF(I,4)-c值;
!
1——垂直于单元的均布力,大小为q,以坐标轴正向为正,c为荷载末端距i节点距离;
!
2——非节点集中力P,c为荷载距i节点距离;
1
计算结构力学程序
!
3——非节点集中力距M,c为荷载距i节点距离,右手法则判正负;
!
4——三角形荷载,c为荷载距i节点距离,i端为0,距离i端c时力为q;
!
j端为0的三角形,可按叠加法处理。
!
5——沿杆轴向均布力,大小为q,c为荷载末端距i节点距离;
!
6——沿杆轴向集中力,大小为q,c为荷载末端距i节点距离;
!
!
从第7到第9类型(支座沉降)数据输法:
PF(I,1)-单元码;PF(I,2)-类型;PF(I,3)-位移大小(含正负),坐
标轴正向位为正,转角按右手法则;PF(I,4)-沉降所在的单元位移分量,i端为1-3,j端为4-6;
!
!
7——沿轴向支座沉降;
!
8——垂直于轴向支座沉降;
!
9——支座转动;
!
10——制造误差,PF(I,1)—制造误差所在单元,PF(I,2)-类型;PF(I,3)-误差大小(含正负),正负取决于消除
!
误差时端点的运动方向,PF(I,4)—误差所在坐标号;
!
11——温度荷载,PF(I,1)—荷载所在单元,数据形式为:
ElementNo.1,如2单元上有温度荷载,则PF(I,1)=2.1;
!
PF(I,2)—温度变化值t1,PF(I,3)—温度变化值t2,PF(I,4)—材料线膨胀系数;
!
!
TK(NN)——采用一维存储结构刚度矩阵,上半带元素(每列第一个非零元素到对角元);
!
KD——主元地址数组,表示结构刚度矩阵的主元在TK中的序号,KD中最后一个数是TK中元素的总个数;
!
JI——结构刚度矩阵上半带的非对角元素在TK中的地址,JI=KD(J)-J+I;
!
JN(NJ,3)——结点位移分量编号数组,用于存放结点三个位移的位移分量号码,
!
JN(I,1),JN(I,2),JN(I,2)-分别为结点I的U,V,CETA分量的位移分量(坐标)号码;
!
!
P(N)——节点荷载列阵;在回代求位移时存放位移量;
!
F(N)——求得的杆端力列阵;
!
FO(6)——等效节点荷载列阵;
!
!
!
!
!
****************************************************************************************
!
!
!
**********************平面结构分析源程序内容**************************************
!
!
!
****************************************************************************************
PROGRAMPFF
DIMENSIONX(50),Y(50),JE(50,2),MT(30),AI(10,2),JS(20,4),PJ(50,3),PF(50,4),JN(50,3),
&KD(150),TK(1000),P(150),F(6),H(50)
DOUBLEPRECISIONTK,P,F
CHARACTER*200TL
OPEN(1,FILE='INDAT.DAT',STATUS='OLD')
OPEN(2,FILE='OUTDAT.DAT',STATUS='NEW')
READ(1,70)TL
READ(1,70)TL
READ(1,*)NE,NJ,NS,NAI,NL,E
2
计算结构力学程序
WRITE(2,10)NE,NJ,NS,NAI,NL,E
10FORMAT(5X,'PLANEFRAMESTRUCTUREANALYSIS'/5X,'**********'//2X,'CONTROLPARAMETERS
&OFSTRUCTURE'/5X,'NE=',I2,8X,'NJ=',I2,8X,'NS=',I2,8X,'NAI=',I2,/5X,'NL=',I2,8X,'E=',E12.4)
CALLINPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H)!
读入数据文件
CALLDJN(NJ,NS,JS,JN,N)!
计算结构自由度数N,形成结点位移分量数组JN
CALLADE(NE,NJ,N,JE,JN,KD,NN)!
形成主元地址数组KD(N)
CALLSSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK)!
形成总刚,一维存储数组TK(NN)
CALLUTDU3(TK,NN,KD,N)!
对总刚进行UTDU分解,以用于解方程组
DO20LC=1,NL!
对工况循环
READ(1,70)TL
READ(1,70)TL
READ(1,70)TL
READ(1,*)NP,NF!
读入工况信息
WRITE(2,30)LC,NP,NF
30FORMAT(/2X,'LOADDATA'/10X,'LOADCASE=',I3/10X,'NP=',I3,8X,'NF=',I3)
CALLNLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H)!
形成总荷载列阵P(N)
CALLBACK3(TK,NN,P,N,KD,JN,NJ)!
解方程组并输出结点位移,存放在数组P(N)中
WRITE(2,40)
40FORMAT(//4X,'MEMBER-ENDFORCESOFELEMENTS'/4X,'ELEMENT',13X,'N',17X,'V',17X,'M')
DO60M=1,NE
CALLMQN(M,NE,NJ,NAI,N,NF,E,X,Y,JE,MT,AI,JN,PF,P,F,H)!
计算单元杆端力,存放在数组F(6)中
WRITE(2,50)M,(F(I),I=1,6)!
输出杆端力
50FORMAT(/1X,I10,3X,'N1=',D12.4,3X,'V1=',D12.4,3X,'M1=',D12.4/14X,'N2=',
&D12.4,3X,'V2=',D12.4,3X,'M2=',D12.4)60CONTINUE
20CONTINUE
70FORMAT(A)
CLOSE
(1)
CLOSE
(2)
END
SUBROUTINEINPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H)!
读入数据文件
DIMENSIONX(NJ),Y(NJ),JE(NE,2),MT(NE),AI(NAI,2),JS(NS,4),H(NE)
INTEGERNO
READ(1,70)TL
READ(1,70)TL
READ(1,70)TL
READ(1,*)(NO,X(I),Y(I),I=1,NJ)
READ(1,70)TL
READ(1,70)TL
READ(1,70)TL
READ(1,*)(NO,JE(I,1),JE(I,2),MT(I),H(I),I=1,NE)
READ(1,70)TL
READ(1,70)TL
3
计算结构力学程序
READ(1,70)TL
READ(1,*)(NO,(AI(I,J),J=1,2),I=1,NAI)
READ(1,70)TL
READ(1,70)TL
READ(1,70)TL
READ(1,*)((JS(I,J),J=1,4),I=1,NS)
WRITE(2,10)(I,X(I),Y(I),I=1,NJ)
WRITE(2,20)(I,JE(I,1),JE(I,2),MT(I),I=1,NE)
WRITE(2,30)(I,(AI(I,J),J=1,2),I=1,NAI)
WRITE(2,40)((JS(I,J),J=1,4),I=1,NS)10FORMAT(//2X,'COORDINATESOFJOINTS'/6X,'JOINT',11X,'X',11X,'Y'/(6X,I4,5X,2F12.4))
20FORMAT(//2X,'INFORMATIONOFELEMENTS'/6X,'ELEMENT',
&4X,'JOINT-I',4X,'JOINT-J',5X,'TYPE'/(2X,4I10))30FORMAT(/7X,'TYPE',10X,'A',12X,'I'/(8X,I2,5X,2F12.6))
40FORMAT(//2X,'INFORMATIONOFSPECIALJOINTS'/6X,'JOINT',4X,'u',4x,'v',4x,'ceta'/(6X,4I5))
70FORMAT(A)
END
!
Determinenumberofjointdisplacements.
SUBROUTINEDJN(NJ,NS,JS,JN,N)!
计算结构自由度数N,形成数组JN。
JN为结点位移分量编
号数组,用于存放结点三个位移的位移分量号码
DIMENSIONJS(NS,4),JN(NJ,3)
DO10I=1,NJ
DO10J=1,3
10JN(I,J)=0!
对JN置0
DO20I=1,NS
L=JS(I,1)!
取特殊节点的节点码
DO20J=1,3
20JN(L,J)=JS(I,J+1)!
将JS中位移约束信息JS(I,2),JS(I,3),JS(I,4)附给JN(L,1),JN(L,2),JN(L,3)
N=0!
自由度置0
ID=0
DO30I=1,NJ
DO30J=1,3
IF(JN(I,J)-1)40,50,60!
由条件的小于0,等于0,大于0来执行40,50,60语句40N=N+1!
小于0时,JN(I,J)=0,有自由度,自由度累加
JN(I,J)=N
GOTO30
50JN(I,J)=0!
等于0,JN(I,J)=1,没有自由度
GOTO30
60ID=1!
大于0,结点I是从结点,自由度不增加,RETURN30CONTINUE
IF(ID.EQ.0)GOTO80!
判断是否有主从结点
DO70I=1,NS
4
计算结构力学程序
L=JS(I,1)
DO70J=1,3
K=JS(I,J+1)
IF(K.LE.1)GOTO70!
若K>1,则L是从结点,K是主结点编号
JN(L,J)=JN(K,J)!
从结点L的第J个位移分量的编号应取主结点K的第J个位移分量的编号70CONTINUE
80RETURN
END
!
Determineaddressofdiagonalelementsoftotalstiffnessmatrix.
SUBROUTINEADE(NE,NJ,N,JE,JN,KD,NN)!
形成主元地址数组KD(N)
DIMENSIONJE(NE,2),JN(NJ,3),KD(N),JC(6)
DO10I=1,N
10KD(I)=0!
主元地址数组KD(N)先置0
DO20M=1,NE!
对每根杆件循环
CALLEJC(M,NE,NJ,JE,JN,JC)!
形成单元定位向量JC
MIN=N!
单元定位向量中的最小非零分量MIN,赋初值为N
DO30I=1,6!
对定位向量的六个分量循环
J=JC(I)!
取定位向量的分量
IF(J.EQ.0)GOTO30!
排除零分量
IF(J.LT.MIN)MIN=J!
求定位向量的最小非零分量
30CONTINUE
DO20I=1,6!
对定位向量的六个分量循环
J=JC(I)!
取定位向量的分量
IF(J.EQ.0)GOTO20!
排除零分量
NW=J-MIN+1!
计算半带宽
IF(NW.GT.KD(J))KD(J)=NW!
第J列的真正半带宽
20CONTINUE
NN=1!
NN置初值为1,NN为存放总刚度矩阵中非零元素的一维数组TK(NN)的元素个数
DO40J=2,N!
对总刚度矩阵[K]列循环
NN=NN+KD(J)
KD(J)=NN!
累加各列半带宽,形成KD数组
40CONTINUE
END
!
Assemblestructurestiffnessmatrixstoredasavector.
SUBROUTINESSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK)!
形成存储总刚的一维数组TK(NN)
DIMENSIONX(NJ),Y(NJ),JE(NE,2),JN(NJ,3),MT(NE),AI(NAI,2),KD(N),TK(NN),KE(6,6),JC(6)
DOUBLEPRECISIONTK,KE,BL,SI,CO
DO10I=1,NN
10TK(I)=0.0!
数组TK(NN)置0
DO20M=1,NE!
对每根杆件循环
CALLSCL(M,NE,NJ,X,Y,JE,BL,SI,CO)!
求单元常数
CALLESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE)!
形成global坐标下的单元刚度矩阵KE(6,6)
5
计算结构力学程序
CALLEJC(M,NE,NJ,JE,JN,JC)!
形成单元定位向量JC
DO30L=1,6!
对单元刚度矩阵的行循环
I=JC(L)!
取单刚的第L行在总刚中的行码
IF(I.EQ.0)GOTO30
DO40K=1,6!
对单元刚度矩阵的列循环
J=JC(K)!
取单刚的第K列在总刚中的列码
IF(J.LT.I)GOTO40!
排除总刚的下三角元素
JI=KD(J)-J+I!
计算总刚第I行第J列元素在数组TK中的地址
TK(JI)=TK(JI)+KE(L,K)!
叠加单刚元素到总刚中40CONTINUE
30CONTINUE
20CONTINUE
!
WRITE(2,50)TK
!
50FORMAT('TK',3X,2D15.4)
END
SUBROUTINESCL(M,NE,NJ,X,Y,JE,BL,SI,CO)!
计算单元常数
DIMENSIONX(NJ),Y(NJ),JE(NE,2)
DOUBLEPRECISIONBL,SI,CO
I=JE(M,1)!
取M单元起始结点号
J=JE(M,2)!
取M单元终止结点号
DX=X(J)-X(I)
DY=Y(J)-Y(I)
BL=SQRT(DX*DX+DY*DY)!
计算单元杆长
SI=DY/BL!
计算单元转角SIN值
CO=DX/BL!
计算单元转角COS值
END
!
DecomposetotalstiffnessmatrixusingCroutmethod.
SUBROUTINEUTDU3(A,NN,ID,N)!
对总刚进行UTDU分解,用于解方程组,对应实参(TK,NN,KD,N)
DIMENSIONA(NN),ID(N)
DOUBLEPRECISIONA
DO30J=2,N!
对列循环
JJ=ID(J)!
取主对角元素在一维数组A(NN)中的地址码,A(NN)对应TK(NN)
J1=J-1
MJ=J-JJ+ID(J1)+1!
第J列第一个非零元素的行号
IF(MJ.GT.J1)GOTO30!
排除下三角元素
DO20I=MJ,J!
对行循环
II=ID(I)
I1=I-1
MIJ=MJ
IF(I.EQ.J)GOTO5
6
计算结构力学程序
MI=I-II+ID(I1)+1
IF(MIJ.LT.MI)MIJ=MI
5S=0.0
IF(MIJ.GT.I1)GOTO15
DO10K=MIJ,I1
KI=II-I+K
KK=ID(K)
KJ=JJ-J+K
10S=S+A(KI)*A(KK)*A(KJ)
15IF(I.EQ.J)THEN
A(JJ)=A(JJ)-S!
求对角阵D的第J个元素,存放在一维数组A(NN)中
ELSE
JI=JJ-J+I
A(JI)=(A(JI)-S)/A(II)!
求三角阵U的第I行第J列元素,存放在一维数组A(NN)中
ENDIF
20CONTINUE
30CONTINUE
END
!
Formtotaljointloadvector.
SUBROUTINENLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H)!
形成总荷载列阵P(N)
DIMENSIONX(NJ),Y(NJ),JE(NE,2),JN(NJ,3),PJ(NP,3),PF(NF,4),MT(NE),AI(NAI,2),
&P(N),FO(6),FE(6),JC(6),H(NE)
DOUBLEPRECISIONP,FO,FE,BL,SI,CO
INTEGERNO
DO10I=1,N
P(I)=0.0!
总荷载列阵P(N)置0
10CONTINUE
READ(1,70)TL
READ(1,70)TL
READ(1,70)TL
IF(NP.GT.0)THEN
READ(1,*)(NO,(PJ(I,J),J=1,3),I=1,NP)!
读入结点荷载信息
WRITE(2,20)((PJ(I,J),J=1,3),I=1,NP)20FORMAT(/2X,'JOINTLOAD'/6X,'JOINT',8X,'XYM',12X,'LOAD'/(6X,F5.0,6X,F5.0,6X,F12.4))
DO30I=1,NP
J=INT(PJ(I,1))!
取结点荷载所在的结点号
K=INT(PJ(I,2))!
取结点荷载在结点上的方向
L=JN(J,K)!
取结点荷载的位移编号
IF(L.NE.0)P(L)=PJ(I,3)!
把结点荷载叠加到总荷载列阵P(N)中30CONTINUE
ENDIF
READ(1,70)TL
READ(1,70)TL
7
计算结构力学程序
READ(1,70)TL
IF(NF.GT.0)THEN
READ(1,*)(NO,(PF(I,J),J=1,4),I=1,NF)!
读入非结点荷载信息
WRITE(2,40)((PF(I,J),J=1,4),I=1,NF)40FORMAT(/2X,'NON-JOINTLOAD'/6X,'ELEMENT',8X,'TYPE',8X,'LOAD',
&12X,'C'/(6X,F6.0,6X,F6.0,4X,2F12.4))
DO50I=1,NF!
对非结点荷载循环
M=INT(PF(I,1))!
取非结点荷载所在的单元码
CALLSCL(M,NE,NJ,X,Y,JE,BL,SI,CO)!
计算单元常数
CALLEFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)!
将非结点荷载转化为等效结点荷载,存放在FO(6)中
CALLEJC(M,NE,NJ,JE,JN,JC)!
形成单元定位数组
FE
(1)=-FO
(1)*CO+FO
(2)*SI!
将等效结点荷载转化到global坐标下,y轴向下的坐标系。
FE
(2)=-FO
(1)*SI-FO
(2)*CO!
将等效结点荷载转化到global坐标下
FE(3)=-FO(3)!
将等效结点荷载转化到global坐标下
FE(4)=-FO(4)*CO+FO(5)*SI!
将等效结点荷载转化到global坐标下
FE(5)=-FO(4)*SI-FO(5)*CO!
将等效结点荷载转化到global坐标下
FE(6)=-FO(6)!
将等效结点荷载转化到global坐标下
DO60J=1,6
L=JC(J)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 计算 结构 力学 fortran 程序