《电算实习报告》word版.docx

结构力学电算实习报告姓名:* 学号:*摘 要: 现实生活中往往存在一些很复杂的连续梁模型，它们的支撑性质和受力特点都比较复杂，结构力学电算则是将它们有规律的进行简化，最后将其所有已知的性质表达为一组数据，将这些数据输入到我们所编好的程序中，既可以得到我们所要求一些量如：杆端力，各截面内力，某界面的内力影响线等。关键词：连续梁、电算程序、截面内力程序功能：连续梁程序可对连续梁结构进行相关的计算。能够计算在任意形式均布荷载、跨中集中力荷载、跨中力偶荷载、节点力偶荷载组合下连续梁单元的杆端力，也可计算连续梁各截面内力，从而绘制连续梁的弯矩图。影响线程序可以计算单位1集中力荷载作用下的指定截面中内力在单位1载荷按指定距离移动情况下的变化情况，从而提取数据绘制影响线。使用说明：（1）连续梁内力程序在连续梁计算程序（.exe文件）所在文件夹中建立一新的txt文件，并命名（1.txt）。在文件其中依次输入单元数、支承类型、节点荷载个数、非节点荷载个数、指定截面个数、弹性模量、杆长&惯性矩（依次按组输入）、节点荷载大小&对应节点编号（依次按组输入）、非节点荷载值&距该单元左端距离&单元号&荷载类型号（依次按组输入）、所求截面所在单元号&距该单元左端距离（依次按组输入），然后保存该文件。运行连续梁内力程序，按提示输入起始文件即初始数据的文件（要求全名），再输入结果文件名，在相应文件夹下打开结果文件即可查看结果。（2）影响线程序在影响线程序（.exe文件）所在文件夹中建立一新的txt文件，并命名。在其中一次输入单元数、支承类型、要求影响线的数量、弹性模量、单位载荷每次移动的距离、杆长&惯性矩、截面所在单元号&距该单元左端距离，然后保存该文件。运行影响线程序，输入起始文件即初始数据的文件，再输入结果文件名，在相应文件夹下打开结果文件即可查看结果。提取期中的数据，利用excel或其它工具绘制影响线即可。程序注释：支承类型：两端铰支为1，左固右铰为2，左铰右固为3，两端固定为4。非节点荷载类型：均布荷载为1，集中力为2，集中力偶为3；节点荷载只有外力偶。非节点荷载作用距离：指到单元左端的距离，其中均布荷载指均布荷载结束位置到单元左端的距离。荷载正负号：力偶逆时针为正，力向上为正。算例1：求图示结构的弯矩图及指定截面内力。【注释】图中形如 的粗线所示断面为任意断面，程序会求出任意断面的剪力及弯矩根据此连续梁的计算简图，将其性质及受力特征表现为一组数据，并输入运行程序进行计算，如下：输入数据:4 4 2 5 7 13 6 5 10 4 8 6 1210 3 4 4-6 3 1 1-5 3 2 2-4 3 3 14 1 3 1-12 3 4 31 1.5 2 1 2 3 3 1 3 3 4 2.99 4 3.01得出如下数据结果输出结果: *连续梁内力计算* 单元数= 4 支承类型= 4 节点荷载个数= 2 非节点荷载个数= 5 指定截面个数= 7 弹性模量= 1.0000 杆长,惯性矩 GC(NE),GX(NE) 3.000 6.000 5.000 10.000 4.000 8.000 6.000 12.000 节点荷载大小,对应未知数序号 PJ(I,1),PJ(I,2) 10.000 3.000 4.000 4.000 非结点荷载值,距离,单元号,荷载类型号 -6.000 3.000 1.000 1.000 -5.000 3.000 2.000 2.000 -4.000 3.000 3.000 1.000 4.000 1.000 3.000 1.000 -12.000 3.000 4.000 3.000 截面所在单元号，距离 1.000 1.500 2.000 1.000 2.000 3.000 3.000 1.000 3.000 3.000 4.000 2.990 4.000 3.010 ：位 移 ：；： 结点号= 1 .0000 结点号= 2 .0109 结点号= 3 .4815 结点号= 4 .5463 结点号= 5 .0000 .各单元杆端内力. 单元号= 1 左端弯矩= 4.543 右端弯矩= -4.413 单元号= 2 左端弯矩= 4.413 右端弯矩= .296 单元号= 3 左端弯矩= 9.704 右端弯矩= 2.630 单元号= 4 左端弯矩= 1.370 右端弯矩= -.815 ;所求截面内力; 截面 1 弯矩= 2.272 剪力= -.043 截面 2 弯矩= -1.471 剪力= -2.942 截面 3 弯矩= 4.412 剪力= -2.942 截面 4 弯矩= -2.621 剪力= -7.083 截面 5 弯矩= 3.546 剪力= .917 截面 6 弯矩= -7.073 剪力= 1.907 截面 7 弯矩= 4.888 剪力= 1.907= 计算结束 =根据此得出数据，绘制连续梁弯矩图：如下：算例2：求图示结构1、2截面的剪力影响线。输入数据：3 1 2 1 14 4 5 5 6 61 2 2 1输出结果： *连续梁弯矩剪力影响线计算* 单元数= 3 支承类型= 1 要求影响线数量= 2 弹性模量= 1.0000 单位荷载每次移动距离 1.0000 杆长,惯性矩 GC(NE),GX(NE) 4.000 4.000 5.000 5.000 6.000 6.000 截面所在单元号，距离 1.000 2.000 2.000 1.000 单位力所在单元号 1 距离 1.0000 截面 1 弯矩= .375 剪力= .313 截面 2 弯矩= -.187 剪力= -.062 单位力所在单元号 1 距离 2.0000 截面 1 弯矩= .800 剪力= -.400 截面 2 弯矩= -.300 剪力= -.100 单位力所在单元号 1 距离 3.0000 截面 1 弯矩= .325 剪力= -.163 截面 2 弯矩= -.262 剪力= -.087 单位力所在单元号 1 距离 4.0000 截面 1 弯矩= .000 剪力= .000 截面 2 弯矩= .000 剪力= .000 单位力所在单元号 2 距离 1.0000 截面 1 弯矩= -.160 剪力= .080 截面 2 弯矩= .512 剪力= -.832 单位力所在单元号 2 距离 2.0000 截面 1 弯矩= -.200 剪力= .100 截面 2 弯矩= .216 剪力= -.616 单位力所在单元号 2 距离 3.0000 截面 1 弯矩= -.160 剪力= .080 截面 2 弯矩= .064 剪力= -.384 单位力所在单元号 2 距离 4.0000 截面 1 弯矩= -.080 剪力= .040 截面 2 弯矩= .008 剪力= -.168 单位力所在单元号 2 距离 5.0000 截面 1 弯矩= .000 剪力= .000 截面 2 弯矩= .000 剪力= .000 单位力所在单元号 3 距离 1.0000 截面 1 弯矩= .051 剪力= -.025 截面 2 弯矩= .000 剪力= .102 单位力所在单元号 3 距离 2.0000 截面 1 弯矩= .074 剪力= -.037 截面 2 弯矩= .000 剪力= .148 单位力所在单元号 3 距离 3.0000 截面 1 弯矩= .075 剪力= -.038 截面 2 弯矩= .000 剪力= .150 单位力所在单元号 3 距离 4.0000 截面 1 弯矩= .059 剪力= -.030 截面 2 弯矩= .000 剪力= .119 单位力所在单元号 3 距离 5.0000 截面 1 弯矩= .032 剪力= -.016 截面 2 弯矩= .000 剪力= .065 单位力所在单元号 3 距离 6.0000 截面 1 弯矩= .000 剪力= .000 截面 2 弯矩= .000 剪力= .000= 计算结束 =1截面弯矩影响线1截面剪力影响线截面2的弯矩剪力影响线同理，不再做出。附：程序清单（程序中含有注释）（1）连续梁内力程序 PROGRAM LXL DIMENSION GC(20),GX(20),PJ(20,2),PF(10,4),sqjm(20,2),DK(2,2),P(45), 1,F0(2),WY(2),F(2),ZK1(45),ZK2(45),gd(20,2),qdnl(20,2) character*12 dat0,dat1 write(*,*)输入:初始数据文件名 read(*,10)dat0 write(*,*)输入:计算结果数据文件名 read(*,10)dat110 format(A12) open(8,file=dat0,status=old) open(9,file=dat1,status=new) write(9,11)11 format(8x,*连续梁内力计算*) READ(8,*)NE,NZ,NP,NF,njm,E0 WRITE(9,12)NE,NZ,NP,NF,njm,E0 12 FORMAT(1X,单元数=,6x,I3,5x,支承类型=,I3,5x,节点荷载个数=I3,/, 1 1x,非节点荷载个数=,I3,5x,指定截面个数=,i3,5x,弹性模量=,F12., 14) NJ=NE+1 CALL SRSJ(NE,NP,NF,GC,GX,PJ,PF,sqjm,njm) CALL XCP(NJ,NP,NF,NE,P,PJ,PF,F0,GC) CALL JCZK(NE,NJ,E0,ZK1,ZK2,DK,GC,GX) CALL ZCCL(NZ,NJ,ZK1,ZK2,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,gd) call jmnl(gd,sqjm,njm,qdnl,gc,ne,nf,pf)ENDC *c 输入初始数据SRSJ子程序C * SUBROUTINE SRSJ(NE,NP,NF,GC,GX,PJ,PF,sqjm,njm) DIMENSION GC(NE),GX(NE),PJ(NP,2),PF(NF,4),sqjm(njm,2)csqjm:所求截面 njm:截面数目 sqjm(njm,1)=所属单元 sqjm(njm,2)=距左端距离C GC(NE)：杆长 GX(NE)：惯性矩C 结点荷载:PJ(NP,1)=荷载大小 PJ(NP,2)=对应未知数序号 C 非结点荷载 : PF(NE,1)=荷载大小 PF(NE,2)=距离 PF(NE,3)=单元号 PF(NE,4)=荷载类型号 WRITE(9,11)11 FORMAT(1X,杆长,惯性矩 GC(NE),GX(NE) READ(8,*)(GC(I),GX(I),I=1,NE) write(9,12)(GC(I),GX(I),I=1,NE)12 FORMAT(10X,2F8.3) IF(NP.GT.0)THEN WRITE(9,22)22 FORMAT(1X,节点荷载大小,对应未知数序号 PJ(I,1),PJ(I,2) READ(8,*)(PJ(I,1),PJ(I,2),I=1,NP) write(9,23)(PJ(I,1),PJ(I,2),I=1,NP)23 FORMAT(1X,2F8.3) END IF IF(NF.GT.0) THEN WRITE(9,33) 33 FORMAT(1X,非结点荷载值,距离,单元号,荷载类型号) READ(8,*)(PF(I,J),J=1,4),I=1,NF) WRITE(9,34)(PF(I,J),J=1,4),I=1,NF) 34 FORMAT(1X,4F8.3) END IFif(njm.gt.0) thenwrite(9,35)35 format(1x,截面所在单元号，距离) read(8,*)(sqjm(i,j),j=1,2),i=1,njm)write(9,36)(sqjm(i,j),j=1,2),i=1,njm)36 format(1x,2f8.3) end if END C *#*c 计算第NHF个非结点荷截引起的等效结点荷截列阵F0C * SUBROUTINE DJH(NHF,NE,NF,PF,F0,GC) DIMENSION PF(NF,4),F0(2),GC(NE)C 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,30),IDC 均布荷载作用下的等效荷载列阵F0(2)10 F0(1)=G*C*C*(6.0-8.0*C1+3.0*C2)/12.0 F0(2)=-G*C*C*C*(4.0-3.0*C1)/12.0/BL GOTO 200C 集中力作用下的等效荷截列阵P0(2)20 F0(1)=G*C*D*D/BL/BL F0(2)=-G*D*C2goto 200C跨中力偶作用下的等效荷载列阵Q0(2) 30F0(1)=+(6.0*G*c*c/bl/bl-6.0*g*c/bl+(2.0*g*c*bl-3.0*g*c*c)/bl/bl+g)F0(2)=-(2.0*g*c*bl-3.0*g*c*c)/bl/bl200 RETURN ENDC *c 计算第NE0个单元的单元刚度矩阵DK(2,2)C * 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)=0.010 CONTINUE15 CONTINUEC 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 ENDc *C 集成总体刚度矩阵,ZKl,ZK2分别存主对角元素和付对角元素c * 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 ENDC *C 形成总荷裁矩阵C * 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 ENDc *c 进行支承条件处理c * SUBROUTINE ZCCL(NZ,NJ,ZK1,ZK2,P) DIMENSION ZK1(NJ),ZK2(NJ),P(NJ) GOTO(10,20,30,20),NZ10 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 ENDc *c 解方程求节点位移Pc * 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 ENDc *c 输出位移c * SUBROUTINE SCWY(NJ,P) DIMENSION P(NJ) WRITE(9,10)10 FORMAT(1X,：位 移 ：；：) DO 100 I=1,NJ WRITE(9,20)I,P(I)20 FORMAT(1X,结点号=,I2,5X,F10.4)100 CONTINUE ENDc *c 计算单元杆端内力c * SUBROUTINE GDL(NE,NJ,NF,E0,P,WY,F,F0,DK,PF,GC,GX,gd) DIMENSION P(NJ),WY(2),F(2),F0(2),DK(2,2),PF(NF,4) DIMENSION GC(NE),GX(NE),gd(ne,2) WRITE(9,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 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)=F(J)-F0(J)40 CONTINUE END IF50 CONTINUE END IF WRITE(9,150) NE0,F(1),F(2)gd(ne0,1)=f(1)!将杆端力数据存入内存 以备调用gd(ne0,2)=f(2)150 FORMAT(1X,单元号=,I2,5X,左端弯矩=,F9.3,2X,右端弯矩=,F9.3)200 CONTINUEendc *c计算任意截面的弯矩和剪力（利用叠加法）c *subroutine jmnl(gd,sqjm,njm,qdnl,gc,ne,nf,pf)dimension gd(ne,2),sqjm(njm,2),gc(ne),qdnl1(nf,njm,2),pf(nf,4),qdn，1l(njm,2)c qdnl(njm,1)=弯矩 qdnl(njm,2)=剪力c 取左侧梁为隔离体，弯矩逆时针为正，剪力向上为正，write(9,20)20format(1x,;所求截面内力;)do njm0=1,njm!按要求的截面顺序循环dyh=sqjm(njm0,1)qdnl(njm0,2)=-(gd(dyh,2)+gd(dyh,1)/gc(dyh)!单元上无荷载情况qdnl(njm0,1)=-(gd(dyh,1)+qdnl(njm0,2)*sqjm(njm0,2)do nf0=1,nfc=pf(nf0,2)d=gc(dyh)-cif(abs(pf(nf0,3)-dyh)=sqjm(njm0,2) thenif(pf(nf0,4)=1)then !单元上有均布荷载且荷载结束位置在所求截面左侧qdnl1(nf0,njm0,2)=pf(nf0,1)*c*(2.0*gc(dyh)-c)/2.0/gc(dyh)-pf(nf0,1,1)*sqjm(njm0,2)qdnl1(nf0,njm0,1)=-(pf(nf0,1)*c*(2.0*gc(dyh)-c)/2.0/gc(dyh)*sqjm(n,1jm0,2)-pf(nf0,1)*sqjm(njm0,2)*sqjm(njm0,2)/2.0)else if(pf(nf0,4)=2)then !单元上有集中力荷载且位置在所求截面左侧qdnl1(nf0,njm0,2)=pf(nf0,1)*d/gc(dyh)qdnl1(nf0,njm0,1)=-qdnl1(nf0,njm0,2)*sqjm(njm0,2)else if(pf(nf0,4)=3)then !单元上有集中力偶且位置在所求截面左侧qdnl1(nf0,njm0,2)=-pf(nf0,1)/gc(dyh)qdnl1(nf0,njm0,1)=-qdnl1(nf0,njm0,2)*sqjm(njm0,2)end if else if(abs(pf(nf0,3)-dyh)1.0e-6).and.(pf(nf0,2)sqjm(njm0,2) &thenif(pf(nf0,4)=1)then !单元上有均布荷载且荷载结束位置在所求截面右侧qdnl1(nf0,njm0,2)=-pf(nf0,1)*c*c/2.0/gc(dyh)qdnl1(nf0,njm0,1)=-pf(nf0,1)*c*c/2.0/gc(dyh)*(gc(dyh)-sqjm(njm0,2),1)else if(pf(nf0,4)=2)then !单元上有集中力荷载且位置在所求截面左侧qdnl1(nf0,njm0,2)=-pf(nf0,1)*c/gc(dyh)qdnl1(nf0,njm0,1)=-pf(nf0,1)*d*c/gc(dyh)-qdnl1(nf0,njm0,2)*(sqjm(n,1jm0,2)-c)else if(pf(nf0,4)=3)then !单元上有集中力偶且位置在所求截面左侧qdnl1(nf0,njm0,2)=-pf(nf0,1)/gc(dyh)qdnl1(nf0,njm0,1)=pf(nf0,1)/gc(dyh)*(sqjm(njm0,2)-gc(dyh)end ifelseqdnl1(nf0,njm0,2)=0 !其他情况下 需叠加的内力为0qdnl1(nf0,njm0,1)=0end ifqdnl(njm0,1)=qdnl(njm0,1)+qdnl1(nf0,njm0,1)!叠加求内力qdnl(njm0,2)=qdnl(njm0,2)+qdnl1(nf0,njm0,2)end do write(9,30) njm0,qdnl(njm0,1),qdnl(njm0,2)30format(1x,截面,i3,2x,弯矩=,f9.3,2x,剪力=,f9.3)end do WRITE(9,11)11 FORMAT(= 计算结束 =)end（2）影响线程序 PROGRAM LXLyxx DIMENSION GC(20),GX(20),PF(1,2),sqjm(20,2),DK(2,2),P(45),F0(2), WY,1(2),F(2),ZK1(45),ZK2(45),gd(20,2),qdnl(20,2),fs(20),fs0(20) character*12 dat0,dat1 write(*,*)输入:初始数据文件名 read(*,10)dat0 write(*,*)输入:计算结果数据文件名 read(*,10)dat110 format(A12) open(8,file=dat0,status=old) open(9,file=dat1,status=new) write(9,11)11 format(8x,*连续梁弯矩剪力影响线计算*) READ(8,*)NE,NZ,njm,E0,hz WRITE(9,12)NE,NZ,njm,E0,hz 12 FORMAT(1X,单元数=,6x,I3,5x,支承类型=,I3,5x,要求影响线数量=, 1i3,5x,弹性模量=,F12.4,6x,单位荷载每次移动距离,f12.4,6x) NJ=NE+1 CALL SRSJ(NE,GC,GX,sqjm,njm) do ne0=1,ne!让单位载荷按指定距离移动的循环fs(ne0)=gc(ne0)/hzfs0(ne0)=int(fs(ne0)+0.1)end done0=0do ne0=1,nepf(1,2)=ne0do fs00=1,fs0(ne0)pf(1,1)=hz*fs00write(9,14) ne0,pf(1,1)14format(/8x,单位力所在单元号,i3,12x,距离,f8.4) CALL XCP(NJ,NE,P,PF,F0,GC) CALL JCZK(NE,NJ,E0,ZK1,ZK2,DK,GC,GX) CALL ZCCL(NZ,NJ,ZK1,ZK2,P) CALL FCQJ(NJ,ZK1,ZK2,P) CALL GDL(NE,NJ,E0,P,WY,F,F0,DK,PF,GC,GX,gd) call jmnl(gd,sqjm,njm,qdnl,gc,ne,pf)end doend do WRITE(9,13)13 FORMAT(= 计算结束 =)ENDC *c 输入初始数据SRSJ子程序C * SUBROUTINE SRSJ(NE,GC,GX,sqjm,njm) DIMENSION GC(NE),GX(NE),sqjm(njm,2)csqjm:所求截面 njm:截面数目 sqjm(njm,1)=所属单元 sqjm(njm,2)=距左端距离C GC(NE)：杆长 GX(NE)：惯性矩C 非结点荷载单位1 : PF(1,1)=距离 PF(1,2)=单元号 WRITE(9,11)11 FORMAT(1X,杆长,惯性矩 GC(NE),GX(NE) READ(8,*)(GC(I),GX(I),I=1,NE) write(9,12)(GC(I),GX(I),I=1,NE)12 FORMAT(10X,2F8.3)if(njm.gt.0) thenwrite(9,35)35 format(1x,截面所在单元号，距离) read(8,*)(sqjm(i,j),j=1,2),i=1,njm)write(9,36)(sqjm(i,j),j=1,2),i=1,njm)36 format(1x,2f8.3) end if END C *#*c 计算单位1力引起的等效结点荷截列阵F0C * SUBROUTINE DJH(NE,PF,F0,GC) DIMENSION PF(1,2),F0(2),GC(NE)C C作用距离,NT-单元号 C=PF(1,1) NT=INT(PF(1,2)+0.1) BL=GC(NT) D=BL-C C1=C/BL C2=C1*C1C 集中力作用下的等效荷截列阵P0(2)20 F0(1)=-1.0*C*D*D/BL/BL F0(2)=1.0*D*C2 ENDC *c 计算第NE0个单元的单元刚度矩阵DK(2,2)C * 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)=0.010 CONTINUE15 CONTINUEC 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 ENDc *C 集成总体刚度矩阵,ZKl,ZK2分别存主对角元素和付对角元素c * 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 ENDC *C 形成总荷裁矩阵C * SUBROUTINE XCP(NJ,NE,P,PF,F0,GC) DIMENSION P(NJ),PF(1,2),F0(2),GC(NE) DO 10 I=1,NJ P(I)=0.010 CONTINUE CALL DJH(NE,PF,F0,GC) ND=INT(PF(1,2)+0.1) P(ND)=P(ND)+F0(1) P(ND+1)=P(ND+1)+F0(2) ENDc *c 进行支承条件处理c * SUBROUTINE ZCCL(NZ,NJ,ZK1,ZK2,P) DIMENSION ZK1(NJ),ZK2(NJ),P(NJ) GOTO(10,20,30,20),NZ10 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 ENDc *c 解方程求节点位移Pc * 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 ENDc *c 输出位移c *c 计算单元杆端内力c * SUBROUTINE GDL(NE,NJ,E0,P,WY,F,F0,DK,PF,GC,GX,gd) DIMENSION P(NJ),WY(2),F(2),F0(2),DK(2,2),PF(1,2) DIMENSION GC(NE),GX(NE),gd(ne,2) 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 CONTINUE35 CONTINUE IF(INT(PF(1,2)+0.1).EQ.NE0) THEN CALL DJH(NE,PF,F0,GC) DO 40 J=1,2 F(J)=F(J)-F0(J)40 CONTINUE END IF50 CONTINUEgd(ne0,1)=f(1)gd(ne0,2)=f(2)200 CONTINUEendc *c计算任意截面的弯矩和剪力c *subroutine jmnl(gd,sqjm,njm,qdnl,gc,ne,pf)dimension gd(ne,2),sqjm(njm,2),gc(ne),qdnl1(1,njm,2),pf(1,2),qdnl(,1njm,2)c qdnl(njm,1)=弯矩 qdnl(njm,2)=剪力c 取左侧梁为隔离体，弯矩逆时针为正，剪力向上为正，do njm0=1,njmdyh=sqjm(njm0,1)qdnl(njm0,2)=-(gd(dyh,2)+gd(dyh,1)/gc(dyh)qdnl(njm0,1)=-(gd(dyh,1)+qdnl(njm0,2)*sqjm(njm0,2)c=pf(1,1)d=gc(dyh)-cif(abs(pf(1,2)-dyh)=sqjm(njm0,2) thenqdnl1(1,njm0,2)=-1.0*d/gc(dyh)qdnl1(1,njm0,1)=-qdnl1(1,njm0,2)*sqjm(njm0,2) else if(abs(pf(1,2)-dyh)1.0e-6).and.(p