结构力学 电算实习

一．完成任务：
1.增加荷载类型（原来程序中只有集中力和均布力）。

在新程序中增加集中力偶和均布力偶两种工况
2.计算指定截面的弯矩和剪力。

3.计算指定截面的弯矩影响线。

4.计算指定截面的剪力影响线。

5.计算指定支座的反力影响线。

二．程序框图：
三程序功能
原始程序利用fortran90编写，程序中变量服从I-N规则。

程序计算仅仅针对于连续梁结构体系，两端支撑类型为刚接或者铰接（梁中间支座均为半铰）。

且默认各单元的抗弯刚度不变，抗拉刚度无穷大，不计剪切变形，弹性模量均相等。

程序中的单元数、节点荷载数、非节点荷载数，可以根据需要通过改变程序的数组大小来实现。

原程序可以求解梁杆端弯矩、节点位移的功能。

可以在梁节点作用弯矩荷载，非节点处仅能作用集中力和均布力荷载。

我在源程序的基础上增加了另外的荷载类型——在非节点出作用中弯矩荷载。

增加了求解杆端剪力的功能。

增加了求解任意指定截面弯矩与剪力的功能。

增加了绘制指定截面弯矩影响线、剪力影响线与指定支座反力影响线的功能。

五．源程序修改部分的说明以及源程序码
! 连续梁静力计算程序
PROGRAM LXL
DIMENSION
GC(20),GX(20),PJ(20,2),PF(10,4),DK(2,2),P(45),F0(2),WY(2),F(2,100),ZK1(45),ZK2(45),Q(2),FJ(2 ,100),FL(100,2),FX(45),WX(45)
WRITE(*,11)
11 FORMAT(1X,'单元数,支承类型,节点荷载,非节点荷载,弹性模量')
open(6,file='初始数据.txt',status='old')
READ(6,*)NE,NZ,NP,NF,E0
NJ=NE+1 !!!节点数比单元数多1
! 输入初始数据
CALL SRSJ(NE,NP,NF,GC,GX,PJ,PF,FL)
! 形成总荷载矩阵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,Q,FJ,WX,FX,FL)
! 求影响线
CALL YXX(NE,NJ,NF,E0,P,WY,F,F0,DK,PF,GC,GX,Q,FJ,WX,FX,FL,ZK1,ZK2,NZ,NP)
close(6)
! 计算结束
END
! ************************************
! 输入初始数据SRSJ子程序
! ************************************
SUBROUTINE SRSJ(NE,NP,NF,GC,GX,PJ,PF,FL)
DIMENSiON GC(NE),GX(NE),PJ(NP,2),PF(NF,4),FL(100,2)
! 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(6,*)(GC(I),GX(I),I=1,NE)
IF(NP.GT.0) THEN
WRITE(*,22)
22 FORMAT(1X,'荷载大小，对应未知数序号 PJ(I,1)，PJ(I,2)')
READ(6,*)(PJ(I,1),PJ(I,2),I=1,NP)
END IF
IF(NF.GT.0) THEN
WRITE(*,33)
33 FORMAT(1X,'非结点荷载值，距离，单元号，荷载类型号')
READ(6,*)((PF(I,J),J=1,4),I=1,NF)
END IF
! 所求截面信息:FL(?,1)=所求截面的单元号 FL(?,2)=所求截面距杆端距离
WRITE(*,44)
44 FORMAT(1X,'所求截面的单元号，距杆端距离，输入完所求截面信息后，输入0,0结束')
I=1
MO=1
DO WHILE(MO>0)
READ(6,*)FL(I,1),FL(I,2)
MO=FL(I,1) !输入完所求截面信息后，输入0,0结束循环。

I=I+1
END DO
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) !!!!q or F
C=PF(NHF,2) !!!!a
NT=INT(PF(NHF,3)+0.1)
ID=INT(PF(NHF,4)+0.1)
BL=GC(NT) !!!!L
D=BL-C !!!!b
C1=C/BL !!!!a/L
C2=D/BL !!!!b/L
GOTO(10,20,30,40), ID
! 均布荷载作用下的等效荷载列阵F0(2)
10 F0(1)=G*C*C*(6-8*C1+3*C1*C1)/12.0
F0(2)=-G*C*C*C*(4-3*C1)/12.0/BL
GOTO 200
! 集中力作用下的等效荷截列阵F0(2)
20 F0(1)=G*C*C2*C2
F0(2)=-G*D*C1*C1
GOTO 200
! 集中力偶作用下的等效荷截列阵F0(2)
30 F0(1)=-G*C2*(2-3*C2)
F0(2)=-G*C1*(2-3*C1)
GOTO 200
!C 4：均布力偶作用下的等效荷截列阵F0(2)
40 F0(1)=G*C*C2*C2
F0(2)=-G*D*C1*C1
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
! ********************************************************* ! 集成总体刚度矩阵，ZK1，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 !!!!1:两端简支，2：左端固定，右端简支，3：左端简支定，右端固定，4：两端固定
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,Q,FJ,WX,FX,FL)
DIMENSION
P(NJ),WY(2),F(2,NE),F0(2),DK(2,2),PF(NF,4),Q(2),FJ(2,NE),WX(45),FX(45),FL(100,2) 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,NE0)=0.0
DO 30 J=1,2
F(I,NE0)=F(I,NE0)+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,NE0)=F(J,NE0)-F0(J)
40 CONTINUE
IF(INT(PF(I,4)+0.1).EQ.1) THEN !将单元等价为简支梁，求出单元在非节点荷载下的杆端剪力。

记为：左端剪力Q(1),右端剪力Q(2)
G=PF(I,1)
C=PF(I,2)
BL=GC(NE0)
D=BL-C
Q(1)=Q(1)-G*C*(BL-C/2.0)/BL
Q(2)=Q(2)+(G*C*C/2.0)/BL
END IF
IF(INT(PF(I,4)+0.1).EQ.2) THEN
G=PF(I,1)
C=PF(I,2)
BL=GC(NE0)
D=BL-C
Q(1)=Q(1)-G*(BL-C)/BL
Q(2)=Q(2)+G*C/BL
END IF
IF(INT(PF(I,4)+0.1).EQ.3) THEN
G=PF(I,1)
C=PF(I,2)
BL=GC(NE0)
D=BL-C
Q(1)=Q(1)+G/BL
Q(2)=Q(2)+G/BL
END IF
50 CONTINUE
END IF
FJ(1,NE0)=(F(1,NE0)+F(2,NE0))/(GC(NE0))+Q(1) !单元在杆端弯矩作用下的杆端剪力加上单元在非节点荷载作用下的杆端剪力即为杆端剪力。

FJ(2,NE0)=(F(2,NE0)+F(1,NE0))/(GC(NE0))+Q(2) !FJ(1,?)=?单元左杆端剪力 FJ(2,?)=?单元右杆端剪力
Q(1)=0
Q(2)=0
WRITE(*,150) NE0,F(1,NE0),F(2,NE0)
150 FORMAT(1X,'单元号=',I2,5X,'左端弯矩=',F9.3,2X,'右端弯矩=',F9.3)
WRITE(*,160) NE0,FJ(1,NE0),FJ(2,NE0)
160 FORMAT(1X,'单元号=',I2,5X,'左端剪力=',F9.3,2X,'右端剪力=',F9.3)
200 CONTINUE
I=1
MO=FL(I,1)
DO WHILE(MO>0)
FX(I)=FJ(1,INT(FL(I,1)+0.1)) !杆端力作用下某截面的内力：FX=截面剪力 WX=截面弯矩
WX(I)=F(1,INT(FL(I,1)+0.1))-FX(I)*FL(I,2)
DO 65 L=1,NF
IF(INT(PF(L,3)+0.1).EQ.INT(FL(I,1)+0.1))THEN !将杆端力作用下某截面的内力与非节点荷载作用下某截面的内力进行叠加，求得某截面内力。

IF(PF(L,2)<=FL(I,2))THEN
IF(INT(PF(L,4)+0.1).EQ.1) THEN
FX(I)=FX(I)+PF(L,1)*PF(L,2)
WX(I)=WX(I)-(FL(I,2)-PF(L,2)+PF(L,2)/2.0)*PF(L,1)*PF(L,2)
END IF
IF(INT(PF(L,4)+0.1).EQ.2) THEN
FX(I)=FX(I)+PF(L,1)
WX(I)=WX(I)-PF(L,1)*(FL(I,2)-PF(L,2))
END IF
IF(INT(PF(L,4)+0.1).EQ.3) THEN
WX(I)=WX(I)+PF(L,1)
END IF
ELSE
IF(INT(PF(L,4)+0.1).EQ.1) THEN
FX(I)=FX(I)+PF(L,1)*FL(I,2)
WX(I)=WX(I)-FL(I,2)/2.0*PF(L,1)*FL(I,2)
END IF
END IF
END IF
65 CONTINUE
WRITE(*,*)'单元号=',INT(FL(I,1)+0.1),'右截面剪力=',FX(I),'右截面弯矩=',WX(I)
MO=FL(I,1)
END DO
END
! *******************************
! 影响线
! *******************************
SUBROUTINE YXX(NE,NJ,NF,E0,P,WY,F,F0,DK,PF,GC,GX,Q,FJ,WX,FX,FL,ZK1,ZK2,NZ,NP)
DIMENSION
P(NJ),WY(2),F(2,NE),F0(2),DK(2,2),PF(NF,4),Q(2),FJ(2,NE),WX(45),FX(45),FL(100,2)
DIMENSION GC(NE),GX(NE),ZK1(NJ),ZK2(NJ)
WRITE(*,55)
55 FORMAT(1X,'求影响线的截面所在的单元，求影响线的截面距杆端的距离，若不求影响线请输入0,0')
READ(6,*)X1,X2
! X1=求影响线的截面所在的单元 X2=求影响线的截面距杆端的距离若不求影响线输入0,0。

IF(X1.GT.0) THEN
NP=0 !结构的节点荷载归零
NF=1 !在结构上作用大小为1的非节点集中力荷载
PJ=0
S=0.01 ! S=步长（单位集中力在结构上移动的步长）
FL=0
FL(1,1)=X1
FL(1,2)=X2
PF=0
PF(1,4)=2
PF(1,1)=-1 !当单位集中力在结构上以S步长移动时，求解单元X1上的距离杆端X2的截面内力
DO U=1,NE
DO V=0,GC(U),S
PF(1,2)=V
PF(1,3)=U
DO 10 I=1,NJ
P(I)=0.0
10 CONTINUE
G=PF(1,1)
C=PF(1,2)
NT=INT(PF(1,3)+0.1)
BL=GC(NT)
D=BL-C
C1=C/BL
C2=C1*C1
F0(1)=G*C*D*D/BL/BL
F0(2)=-G*D*C2
P(NT)=P(NT)+F0(1)
P(NT+1)=P(NT+1)+F0(2)
CALL JCZK(NE,NJ,E0,ZK1,ZK2,DK,GC,GX)
CALL ZCCL(NZ,NJ,ZK1,ZK2,P)
CALL FCQJ(NJ,ZK1,ZK2,P)
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,NE0)=0.0
DO 30 J=1,2
F(I,NE0)=F(I,NE0)+DK(I,J)*WY(J)
30 CONTINUE
35 CONTINUE
Q(1)=0
Q(2)=0
IF(INT(PF(1,3)+0.1).EQ.NE0) THEN
G=PF(1,1)
C=PF(1,2)
NT=INT(PF(1,3)+0.1)
BL=GC(NT)
D=BL-C
C1=C/BL
C2=C1*C1
F0(1)=G*C*D*D/BL/BL
F0(2)=-G*D*C2
DO 40 J=1,2
F(J,NE0)=F(J,NE0)-F0(J)
40 CONTINUE
Q(1)=-G*(BL-C)/BL
Q(2)=G*C/BL
END IF
FJ(1,NE0)=(F(1,NE0)+F(2,NE0))/(GC(NE0))+Q(1) FJ(2,NE0)=(F(2,NE0)+F(1,NE0))/(GC(NE0))+Q(2) 200 CONTINUE
FX(1)=FJ(1,INT(FL(1,1)+0.1))
WX(1)=F(1,INT(FL(1,1)+0.1))-FX(1)*FL(1,2)
IF(INT(PF(1,3)+0.1).EQ.INT(FL(1,1)+0.1))THEN IF(PF(1,2)<=FL(1,2))THEN
FX(1)=FX(1)+PF(1,1)
WX(1)=WX(1)-PF(1,1)*(FL(1,2)-PF(1,2))
END IF
END IF
open(9,file='影响线.txt',status='replace')
! 分别输出：影响线横坐标，剪力影响线的值，弯矩影响线的值
write(9,'(1X,F9.3,F9.3,F9.3)')(V+U-1),FX(1),-WX(1)
open(8,file='支座影响线.txt')
! 分别输出：影响线横坐标，左端支座弯矩影响线的值，左端支座剪力影响线的值，右端支座弯矩影响线的值，右端支座剪力影响线的值
write(8,'(1X,F9.3,F9.3,F9.3,F9.3,F9.3)')(V+U-1),F(1,1),FJ(1,1),F(2,NE),-FJ(2,NE)
END DO
END DO
close(9)
close(8)
END IF
WRITE(*,11)
11 FORMAT(1X,'====================== 计算结束 ======================')
END
六．算例：以三跨连续梁为例
两端简支，单元长度为1，惯性矩EI=1，单元1,3的跨中作用大小为1的单位力，单元2作用大小为1的均布力，利用FORTRAN程序进行计算并和结构力学求解器所得结果进行对比。

输入初始数据：
!!!!单元数,支承类型,节点荷载,非节点荷载,弹性模量
3 1 0 3 1
!!!!杆长，惯性矩 GC(NE),GX(NE)
10 1 1 1 1 1
!!!!非结点荷载值，距离，单元号，荷载类型号
12 0.5 1 2
12 1 2 1
12 0.5 3 2
!!!!所求截面的单元号，距杆端距离输入完所求截面信息后，输入0,0结束
1 0.5
0,0
!!!!求影响线的截面所在的单元，求影响线的截面距杆端的距离，若不求影响线请输入0,0
1,0.5
输出结果：
：：：：：：：：：位移：；：：：：：：：：
节点号1 0.4417
节点号2 -0.5208
节点号3 0.5208
节点号4 -0.4417
..............各单元杆端内力................
单元号1 左端弯矩=0.000 右端弯矩=0.625
单元号1 左端剪力=-0.375 右端剪力=0.625
单元号2 左端弯矩=-0.625 右端弯矩=0.825
单元号2 左端剪力=-0.500 右端剪力=0.500
单元号3 左端弯矩=-0.825 右端弯矩=0.000
单元号3 左端剪力=-0.625 右端剪力=0.375
（影响线的计算结果在该文件夹中，因为数据过于庞大，在这里不给予复制）与结构力学求解器所得结果进行对比
弯矩图
剪力图
由此可知，程序计算所得单元内力结果与结构力学求解器所得结果相差无几绘制剪力弯矩影响线（所选截面为单元1跨中截面）
利用程序结算所得数据绘制剪力及弯矩影响线：
剪力影响线：
弯矩影响线：
由此可知，程序计算所得单元内力结果与结构力学求解器所得结果相差无几
绘制支座反力影响线：
左端支座剪力影响线：（弯矩影响线数值为0）
右端支座剪力影响线：（弯矩影响线数值为0）
由于结构力学求解器无法求出支座反力影响线，故不作对比
七．实习体会与结论综述
经过为期一周的电算实习，我通过对连续梁计算程序的改进与应用我更能熟悉利用fortran编译软件解决问题，同时也进一步加深了自己对矩阵位移法的理解，学到了不少的知识。

也学会了一些其他有用的软件例如：EXCEL、求解器等。

这些都对我以后的学习与工作有很大的提高。

实习时间虽短但我获益匪浅，我觉得以后我们应该多组织一些这样的实习，使我们在各方面都能得到很好的提高。