连续梁计算程序

! 连续梁静力计算程序
PROGRAM LXL
DIMENSION GC(20),GX(20),PJ(20,2),PF(10,4),DK(2,2),P(45),F0(2),WY(2),F(2),ZK1(45),ZK2(45)
WRITE(*,11)
11 FORMAT(1X,'单元数,支承类型,节点荷载,非节点荷载,弹性模量')
READ(*,*)NE,NZ,NP,NF,E0
NJ=NE+1
! 输入初始数据
CALL SRSJ(NE,NP,NF,GC,GX,PJ,PF)
! 形成总荷载矩阵P
CALL XCP(NJ,NP,NF,NE,P,PJ,PF,F0,GC)
! 集成整体刚度矩阵ZKl，ZK2
CALL JCZK(NE,NJ,E0,ZK1,ZK2,DK,GC,GX)
! 支承条件处理
CALL ZCCL(NZ,NJ,ZK1,ZK2,P)
! 方程求解--求位移P
CALL FCQJ(NJ,ZK1,ZK2,P)
! 输出位移
CALL SCWY(NJ,P)
! 计算杆端内力
CALL GDL(NE,NJ,NF,E0,P,WY,F,F0,DK,PF,GC,GX)
! 计算结束
END
! ************************************
! 输入初始数据SRSJ子程序
! ************************************
SUBROUTINE SRSJ(NE,NP,NF,GC,GX,PJ,PF)
DIMENSiON GC(NE),GX(NE),PJ(NP,2),PF(NF,4)
! GC(NE)：杆长 GX(NE)：惯性矩
! 结点荷载:PJ(NP,1)=荷载大小 PJ(NP,2)=对应未知数序号
! 非结点荷载 : PF(NE,1)=荷载大小 PF(NE,2)=距离 PF(NE,3)=单元号 PF(NE,4)=荷载类型号
WRITE(*,11)
11 FORMAT(1X,'杆长，惯性矩 GC(NE),GX(NE)')
READ(*,*)(GC(I),GX(I),I=1,NE)
IF(NP.GT.0) THEN
WRITE(*,22)
22 FORMAT(1X,'荷载大小，对应未知数序号 PJ(I,1)，PJ(I,2)')
READ(*,*)(PJ(I,1),PJ(I,2),I=1,NP)
END IF
IF(NF.GT.0) THEN
WRITE(*,33)
33 FORMAT(1X,'非结点荷载值，距离，单元号+荷载类型号')
READ(*,*)((PF(I,J),J=1,4),I=1,NF)
END IF
END
! *******************************#***************
! 计算第NHF个非结点荷截引起的等效结点荷截列阵F0
! ***********************************************
SUBROUTINE DJH(NHF,NE,NF,PF,F0,GC)
DIMENSION PF(NF,4),F0(2),GC(NE)
! G—荷载值，C—作用距离，NT--单元号，ID—荷载类型号
G=PF(NHF,1)
C=PF(NHF,2)
NT=INT(PF(NHF,3)+0.1)
ID=INT(PF(NHF,4)+0.1)
BL=GC(NT)
D=BL-C
C1=C/BL
C2=C1*C1
GOTO(10,20), ID
! 均布荷载作用下的等效荷载列阵F0(2)
10 F0(1)=G*C*C*(6.-8.*C1+3.*C2)/12.0
F0(2)=-G*C*C*C*(4.-3.*C1)/12.0/BL
GOTO 200
! 集中力作用下的等效荷截列阵P0(2)
20 F0(1)=G*C*D*D/BL/BL
F0(2)=-G*D*C2
200 RETURN
END
! ******************************************
! 计算第NE0个单元的单元刚度矩阵DK(2，2)
! ******************************************
SUBROUTINE DG(NE0,NE,E0,GC,GX,DK)
DIMENSION GC(NE),GX(NE),DK(2,2)
DO 15 I=1,2
DO 10 J=1,2
DK(I,J)=1.0
10 CONTINUE
15 CONTINUE
! DL—杆长，DI-惯性矩，S--线刚度
DL=GC(NE0)
DI=GX(NE0)
S=E0*DI/DL
DK(1,1)=4.0*S

DK(1,2)=2.0*S
DK(2,1)=2.0*S
DK(2,2)=4.0*S
END
! *********************************************************
! 集成总体刚度矩阵，ZKl，ZK2分别存主对角元素和付对角元素
! *********************************************************
SUBROUTINE JCZK(NE,NJ,E0,ZK1,ZK2,DK,GC,GX)
DIMENSION ZK1(NJ),ZK2(NJ),DK(2,2),GC(NE),GX(NE)
DO 10 I=1,NJ
ZK1(I)=0.0
ZK2(I)=0.0
10 CONTINUE
DO 100 I=1,NE
CALL DG(I,NE,E0,GC,GX,DK)
ZK1(I)=ZK1(I)+DK(1,1)
ZK2(I)=ZK2(I)+DK(1,2)
ZK1(I+1)=ZK1(I+1)+DK(2,2)
100 CONTINUE
END
! *****************************
! 形成总荷裁矩阵
! *****************************
SUBROUTINE XCP(NJ,NP,NF,NE,P,PJ,PF,F0,GC)
DIMENSION P(NJ),PJ(NP,2),PF(NF,4),F0(2),GC(NE)
DO 10 I=1,NJ
P(I)=0.0
10 CONTINUE
IF(NP.GT.0) THEN
DO 20 I=1,NP
J=INT(PJ(I,2)+0.1)
P(J)=PJ(I,1)
20 CONTINUE
END IF
IF(NF.GT.0) THEN
DO 100 NF0=1,NF
CALL DJH(NF0,NE,NF,PF,F0,GC)
ND=INT(PF(NF0,3)+0.1)
P(ND)=P(ND)+F0(1)
P(ND+1)=P(ND+1)+F0(2)
100 CONTINUE
END IF
END
! **********************************
! 进行支承条件处理
! **********************************
SUBROUTINE ZCCL(NZ,NJ,ZK1,ZK2,P)
DIMENSION ZK1(NJ),ZK2(NJ),P(NJ)
GOTO(10,20,30,20),NZ
10 GOTO 100
20 ZK1(1)=1.0
P(1)=0.0
ZK2(1)=0.0
IF(NZ.EQ.4) GOTO 30
GOTO 100
30 ZK1(NJ)=1.0
ZK2(NJ-1)=0.0
P(NJ)=0.0
100 CONTINUE
END
! ******************************
! 解方程求节点位移P
! ******************************
SUBROUTINE FCQJ(NJ,ZK1,ZK2,P)
DIMENSION ZK1(NJ),ZK2(NJ),P(NJ)
DO 10 I=1,NJ-1
ZK1(I+1)=ZK1(I+1)-ZK2(I)*ZK2(I)/ZK1(I)
P(I+1)=P(I+1)-ZK2(I)*P(I)/ZK1(I)
10 CONTINUE
P(NJ)=P(NJ)/ZK1(NJ)
DO 20 I=1,NJ-1
P(NJ-I)=(P(NJ-I)-ZK2(NJ-I)*P(NJ-I+1))/ZK1(NJ-I)
20 CONTINUE
END
! *******************************
! 输出位移
! *******************************
SUBROUTINE SCWY(NJ,P)
DIMENSION P(NJ)
WRITE(*,10)
10 FORMAT(1X,'：：：：：：：：：位 移 ：；：：：：：：：：')
DO 100 I=1,NJ
WRITE(*,20)I,P(I)
20 FORMAT(1X,'结点号=',I2,5X,F10.4)
100 CONTINUE
END
! ******************************
! 计算单元杆端内力
! ******************************
SUBROUTINE GDL(NE,NJ,NF,E0,P,WY,F,F0,DK,PF,GC,GX)
DIMENSION P(NJ),WY(2),F(2),F0(2),DK(2,2),PF(NF,4)
DIMENSION GC(NE),GX(NE)
WRITE(*,2)
2 FORMAT(1X,'.................各单元杆端内力....................')
DO 200 NE0=1,NE
CALL DG (NE0,NE,E0,GC,GX,DK)
WY(1)=P(NE0)
WY(2)=P(NE0+1)
DO 35 I=1,2
F(I)=0.0
DO 30 J=1,2
F(I)=F(I)

+DK(I,J)*WY(J)
30 CONTINUE
35 CONTINUE
IF(NF.GT.0) THEN
DO 50 I=1,NF
IF(INT(PF(I,3)+0.1).EQ.NE0) THEN
CALL DJH(I,NE,NF,PF,F0,GC)
DO 40 J=1,2
F(J)=F(J)-F0(J)
40 CONTINUE
END IF
50 CONTINUE
END IF
WRITE(*,150) NE0,F(1),F(2)
150 FORMAT(1X,'单元号=',I2,5X,'左端弯矩=',F9.3,2X,'右端弯矩',F9.3)
200 CONTINUE
WRITE(*,11)
11 FORMAT(1X,'====================== 计算结束 ======================')
END