结构力学 电算实习.doc

一完成任务：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(MO0) 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.010 CONTINUE15 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.010 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.010 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 10020 ZK1(1)=1.0 P(1)=0.0 ZK2(1)=0.0 IF(NZ.EQ.4) GOTO 30 GOTO 10030 ZK1(NJ)=1.0 ZK2(NJ-1)=0.0 P(NJ)=0.0100 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 CONTINUE35 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 END IF50 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(MO0) 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 IF65 CONTINUE WRITE(*,*)单元号=,INT(FL(I,1)+0.1),右截面剪力=,FX(I),右截面弯矩=,WX(I) I=I+1 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.010 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 CONTINUE35 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的均布