有限元计算结构力学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）
升级会员 