有限元计算结构力学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)——等效节点荷载列阵;

!

!

!!!*****************************************************************

***********************

!!!********************** 平面结构分析源程序内容

**************************************

!!!*****************************************************************

***********************

PROGRAM PFF

DIMENSION

X(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)

DOUBLE PRECISION TK,P,F

CHARACTER *200 TL

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

10 FORMAT(5X,'PLANE FRAME STRUCTURE

ANALYSIS'/5X,'**********'//2X,'CONTROL PARAMETERS

&OF

STRUCTURE'/5X,'NE=',I2,8X,'NJ=',I2,8X,'NS=',I2,8X,'NAI=',I2,/5X,'NL=',I2 ,8X,'E=',E12.4)

CALL INPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H) !读入数据文件

CALL DJN(NJ,NS,JS,JN,N) !计算结构自由度数N，形成结点位移分量数组JN CALL ADE(NE,NJ,N,JE,JN,KD,NN) !形成主元地址数组KD(N)

CALL SSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK) !形成总刚，一维存储数组TK(NN)

CALL UTDU3(TK,NN,KD,N) !对总刚进行UTDU分解，以用于解方程组

DO 20 LC=1,NL !对工况循环

READ(1,70)TL

READ(1,70)TL

READ(1,70)TL

READ(1,*)NP,NF !读入工况信息

WRITE(2,30)LC,NP,NF

30 FORMAT(/2X,'LOAD DATA'/10X,'LOAD

CASE=',I3/10X,'NP=',I3,8X,'NF=',I3)

CALL NLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H) !形成总荷载列阵P(N)

CALL BACK3(TK,NN,P,N,KD,JN,NJ) !解方程组并输出结点位移，存放在数组

P(N)中

WRITE(2,40)

40 FORMAT(//4X,'MEMBER-END FORCES OF

ELEMENTS'/4X,'ELEMENT',13X,'N',17X,'V',17X,'M')

DO 60 M=1,NE

CALL MQN(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) !输出杆端力

50

FORMAT(/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) 60 CONTINUE

20 CONTINUE

70 FORMAT(A)

CLOSE(1)

CLOSE(2)

END

SUBROUTINE INPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H) !读入数据文件

DIMENSION X(NJ),Y(NJ),JE(NE,2),MT(NE),AI(NAI,2),JS(NS,4),H(NE)

INTEGER NO

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) 10 FORMAT(//2X,'COORDINATES OF JOINTS'/6X,'JOINT',11X,'X',11X,'Y'/(6X,I4,5X,2F12.4))

20 FORMAT(//2X,'INFORMATION OF ELEMENTS'/6X,'ELEMENT',

& 4X,'JOINT-I',4X,'JOINT-J',5X,'TYPE'/(2X,4I10)) 30

FORMAT(/7X,'TYPE',10X,'A',12X,'I'/(8X,I2,5X,2F12.6))

40 FORMAT(//2X,'INFORMATION OF SPECIAL

JOINTS'/6X,'JOINT',4X,'u',4x,'v',4x,'ceta'/(6X,4I5))

70 FORMAT(A)

END

! Determine number of joint displacements.

SUBROUTINE DJN(NJ,NS,JS,JN,N) !计算结构自由度数N，形成数组JN。JN为结点位移分量编

号数组，用于存放结点三个位移的位移分量号码

DIMENSION JS(NS,4),JN(NJ,3)

DO 10 I=1,NJ

DO 10 J=1,3

10 JN(I,J)=0 !对JN置0

DO 20 I=1,NS

L=JS(I,1) !取特殊节点的节点码

DO 20 J=1,3

20 JN(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

DO 30 I=1,NJ

DO 30 J=1,3

IF(JN(I,J)-1)40,50,60 !由条件的小于0，等于0，大于0来执行40，50，60语句 40 N=N+1 !小于0时，JN(I,J)=0，有自由度，自由度累加

JN(I,J)=N

GOTO 30

50 JN(I,J)=0 !等于0，JN(I,J)=1，没有自由度

GOTO 30

60 ID=1 !大于0，结点I是从结点,自由度不增加，RETURN 30 CONTINUE

IF(ID.EQ.0)GOTO 80 !判断是否有主从结点

DO 70 I=1,NS

4

计算结构力学程序

L=JS(I,1)

DO 70 J=1,3

K=JS(I,J+1)

IF(K.LE.1)GOTO 70 !若K>1，则L是从结点，K是主结点编号

JN(L,J)=JN(K,J) !从结点L的第J个位移分量的编号应取主结点K的第J个位移分量的编号 70 CONTINUE

80 RETURN

END

! Determine address of diagonal elements of total stiffness matrix.

SUBROUTINE ADE(NE,NJ,N,JE,JN,KD,NN) !形成主元地址数组KD(N)

DIMENSION JE(NE,2),JN(NJ,3),KD(N),JC(6)

DO 10 I=1,N

10 KD(I)=0 !主元地址数组KD(N)先置0

DO 20 M=1,NE !对每根杆件循环

CALL EJC(M,NE,NJ,JE,JN,JC) !形成单元定位向量JC

MIN=N !单元定位向量中的最小非零分量MIN，赋初值为N

DO 30 I=1,6 !对定位向量的六个分量循环

J=JC(I) !取定位向量的分量

IF(J.EQ.0) GOTO 30 !排除零分量

IF(J.LT.MIN) MIN=J !求定位向量的最小非零分量

30 CONTINUE

DO 20 I=1,6 !对定位向量的六个分量循环

J=JC(I) !取定位向量的分量

IF(J.EQ.0) GOTO 20 !排除零分量

NW=J-MIN+1 !计算半带宽

IF(NW.GT.KD(J)) KD(J)=NW !第J列的真正半带宽

20 CONTINUE

NN=1 !NN置初值为1，NN为存放总刚度矩阵中非零元素的一维数组TK(NN)的

元素个数

DO 40 J=2,N !对总刚度矩阵[K]列循环

NN=NN+KD(J)

KD(J)=NN !累加各列半带宽，形成KD数组

40 CONTINUE

END

! Assemble structure stiffness matrix stored as a vector.

SUBROUTINE SSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK) !形成存储总刚

的一维数组TK(NN)

DIMENSION

X(NJ),Y(NJ),JE(NE,2),JN(NJ,3),MT(NE),AI(NAI,2),KD(N),TK(NN),KE(6,6),JC(6) DOUBLE PRECISION TK,KE,BL,SI,CO

DO 10 I=1,NN

10 TK(I)=0.0 !数组TK(NN)置0

DO 20 M=1,NE !对每根杆件循环

CALL SCL(M,NE,NJ,X,Y,JE,BL,SI,CO) !求单元常数

CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成global坐标下的单元刚度矩阵KE(6,6)

5

计算结构力学程序

CALL EJC(M,NE,NJ,JE,JN,JC) !形成单元定位向量JC

DO 30 L=1,6 !对单元刚度矩阵的行循环

I=JC(L) !取单刚的第L行在总刚中的行码

IF(I.EQ.0) GOTO 30

DO 40 K=1,6 !对单元刚度矩阵的列循环

J=JC(K) !取单刚的第K列在总刚中的列码

IF(J.LT.I) GOTO 40 !排除总刚的下三角元素

JI=KD(J)-J+I !计算总刚第I行第J列元素在数组TK中的地址

TK(JI)=TK(JI)+KE(L,K) !叠加单刚元素到总刚中 40 CONTINUE

30 CONTINUE

20 CONTINUE

! WRITE(2,50) TK

!50 FORMAT('TK',3X,2D15.4)

END

SUBROUTINE SCL (M,NE,NJ,X,Y,JE,BL,SI,CO) !计算单元常数

DIMENSION X(NJ),Y(NJ),JE(NE,2)

DOUBLE PRECISION BL,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

! Decompose total stiffness matrix using Crout method.

SUBROUTINE UTDU3(A,NN,ID,N)!对总刚进行UTDU分解，用于解方程组，对应实参(TK,NN,KD,N)

DIMENSION A(NN),ID(N)

DOUBLE PRECISION A

DO 30 J=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) GOTO 30 !排除下三角元素

DO 20 I=MJ,J !对行循环

II=ID(I)

I1=I-1

MIJ=MJ

IF(I.EQ.J) GOTO 5

6

计算结构力学程序

MI=I-II+ID(I1)+1

IF(MIJ.LT.MI) MIJ=MI

5 S=0.0

IF(MIJ.GT.I1) GOTO 15

DO 10 K=MIJ,I1

KI=II-I+K

KK=ID(K)

KJ=JJ-J+K

10 S=S+A(KI)*A(KK)*A(KJ)

15 IF(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)中

END IF

20 CONTINUE

30 CONTINUE

END

! Form total joint load vector.

SUBROUTINE NLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H) !形成总荷载列阵P(N)

DIMENSION

X(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)

DOUBLE PRECISION P,FO,FE,BL,SI,CO

INTEGER NO

DO 10 I=1,N

P(I)=0.0 !总荷载列阵P(N)置0

10 CONTINUE

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) 20 FORMAT(/2X,'JOINT LOAD'/6X,'JOINT',8X,'XYM',12X,'LOAD'/(6X,F5.0,6X,F5.0,6X,F12.4)) DO 30 I=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)中 30 CONTINUE

END IF

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) 40 FORMAT(/2X,'NON-JOINT

LOAD'/6X,'ELEMENT',8X,'TYPE',8X,'LOAD',

&12X,'C'/(6X,F6.0,6X,F6.0,4X,2F12.4))

DO 50 I=1,NF !对非结点荷载循环

M=INT(PF(I,1)) !取非结点荷载所在的单元码

CALL SCL(M,NE,NJ,X,Y,JE,BL,SI,CO) !计算单元常数

CALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)!将非结点荷载转化为等效结点荷载，存放在FO(6)中

CALL EJC(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坐标下

DO 60 J=1,6

L=JC(J)

IF(L.NE.0) P(L)=P(L)+FE(J) !将等效结点荷载叠加到总荷载列阵P(N)中 60 CONTINUE

50 CONTINUE

END IF

70 FORMAT(A)

END

! Solve equations and print joint displacements.

SUBROUTINE BACK3(A,NN,B,N,ID,JN,NJ) !解方程组并输出结点位移，存放在数组B(N)中。这里B(N)

数组对应主程序中的P(N)数组

DIMENSION A(NN),B(N),ID(N),JN(NJ,3),D(3)

DOUBLE PRECISION A,D,B

DO 20 J=2,N !对列循环

JJ=ID(J)

J1=J-1

MJ=J-JJ+ID(J1)+1 !第J列第一个非零元素的行号

IF(MJ.GT.J1) GOTO 20 !若MJ>J-1，则z的第J个值即为B的第J个值

DO 10 I=MJ,J1

JI=JJ-J+I

10 B(J)=B(J)-A(JI)*B(I) !求z的第J个值，存放在B(N)数组。这里B(N)数组对应主程序中的P(N)数组 20 CONTINUE

DO 30 I=1,N

II=ID(I)

B(I)=B(I)/A(II) !求y的第J个值，存放在B(N)数组

30 CONTINUE

8

计算结构力学程序

NW=1 !半带宽置初值为1

DO 35 I=2,N

IW=ID(I)-ID(I-1)

IF(IW.GT.NW) NW=IW !计算最大列半带宽

35 CONTINUE

N1=N-1

DO 50 I=N1,1,-1 !对行循环

MIN=I+NW-1

IF(N.LT.MIN) MIN=N

DO 40 J=I+1,MIN !对列循环

JJ=ID(J)

MJ=J-JJ+ID(J-1)+1 !第J列第一个非零元素的行号

IF(I.LT.MJ) GOTO 40 !判断是否为带内元素

JI=JJ-J+I

B(I)=B(I)-A(JI)*B(J) !求x的第I个值，存放在B(N)数组。这里x的值即为结点位移

40 CONTINUE

50 CONTINUE

WRITE(2,60)

60 FORMAT(//2X,'JOINT

DISPLACEMENTS'/5X,'JOINT',12X,'U',14X,'V',11X,'ceta')

DO 80 I=1,NJ

DO 70 J=1,3

D(J)=0.0

L=JN(I,J) !取位移号

IF(L.NE.0) D(J)=B(L) !给位移赋值

70 CONTINUE

WRITE(2,75) I,D(1),D(2),D(3) !输出第I号结点的位移值 75

FORMAT(2X,I6,4X,3D15.6) 80 CONTINUE

END

SUBROUTINE MQN(M,NE,NJ,NAI,N,NF,E,X,Y,JE,MT,AI,JN,PF,P,F,H)!计算单元杆端力，存放在数组F(6)中

DIMENSION

X(NJ),Y(NJ),JE(NE,2),AI(NAI,2),P(N),JN(NJ,3),PF(NF,4),JC(6),FO(6),D(6),F D(6),

&F(6),KE(6,6),MT(NE),H(NE)

DOUBLE PRECISION P,FO,D,FD,F,KE,BL,SI,CO

CALL SCL(M,NE,NJ,X,Y,JE,BL,SI,CO) !求单元常数

CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成global坐标下的单元刚度矩阵KE(6,6)

CALL EJC(M,NE,NJ,JE,JN,JC) !形成单元定位向量JC

DO 10 I=1,6

L=JC(I) !取位移编号

D(I)=0.0

IF(L.NE.0) D(I)=P(L) !取global坐标下的杆端位移 10 CONTINUE

DO 20 I=1,6

FD(I)=0.0

9

计算结构力学程序

DO 30 J=1,6

FD(I)=FD(I)+KE(I,J)*D(J) !计算杆端力

30 CONTINUE

20 CONTINUE

F(1)=FD(1)*CO+FD(2)*SI !将杆端力转化到local坐标下

F(2)=-FD(1)*SI+FD(2)*CO !将杆端力转化到local坐标下

F(3)=FD(3) !将杆端力转化到local坐标下

F(4)=FD(4)*CO+FD(5)*SI !将杆端力转化到local坐标下

F(5)=-FD(4)*SI+FD(5)*CO !将杆端力转化到local坐标下

F(6)=FD(6) !将杆端力转化到local坐标下

IF(NF.GT.0)THEN

DO 40 I=1,NF

L=INT(PF(I,1))

IF(M.EQ.L)THEN

CALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)

DO 50 J=1,6

F(J)=F(J)+FO(J) !如果有非结点荷载，杆端力中要加入等效结点荷载 50 CONTINUE

ENDIF

40 CONTINUE