空间刚架计算程序及说明.doc

一、 设计意图及总体框架：大多数空间刚架都属于超静定结构，因此难于用一般的结构力学方法进行计算，如用矩阵位移法则由于阶数太大难于手算，但如果结合有限元方法编写程序进行计算，则会有很大的普遍性，且方便快捷。本程序能计算只承受任意的非节点集中力、均布荷载和集中力偶矩的所有刚架（包括平面刚架、平面板架）的节点内力和位移，但不适用有支座沉降和温度变化情况，计算出节点力之后，可以很方便地画出结构内力图，从而整个问题都得到求解。程序设计的整个过程可以参阅相关书籍，空间刚架的程序设计和平面刚架的程序设计的最大区别在于空间刚架计算的坐标变换矩阵较复杂和由于维数和阶数的增加所带来的子程序之间衔接的复杂性，除此之外，一些平面刚架程序的子函数可以直接用到空间刚架程序中来。本程序一共包括10子函数：局部坐标下单元刚度计算子程序KEP()、坐标变换矩阵计算子程序CR()、固端反力计算子程序FIXF()、单元定位向量计算子程序（）、指示数组计算子程序LD()、约束处理子程序FCC()、解线性方程组子程序DECOM()、矩阵赋零子程序ZERO()、矩阵转置子程序TRAN()和矩阵相乘子程序DOT()。具体算法见源程序清单，输入输出都用文件进行。注：1,为了方便求坐标变换矩阵，可以在建立局部坐标系时使y/XOY,而这一点我们总是可以做到的 2,在上面的坐标系下，可以将刚架单元分为竖直单元和一般单元，分别求坐标转换矩阵，具体做法见程序。 3,单元刚度矩阵可以用平面刚架和平面板架的结果叠加。二、 变量及符号：NN节点总数， NM材料类型数 NA截面类型数NE单元总数， NC受约束自由度数 NF节点自由度数NP2非节点荷载数 ND单元节点数 NFD单元自由度数N结构自由度总数 NT一维数组A的容量 X(NN)X坐标数组， Y(NN)Y坐标数组 Z(NN)Z坐标数组MEA(NE,4)单元信息数组，其中: MEA(I,1)I单元i节点编号 MEA(I,2)I单元j节点编号 MEA(I,3)I单元截面编号 MEA(I,4)I单元材料编号AEU(NM,3)材料信息数组，其中: AEU(NM1,1)弹性模量EO(E) AEU(NM1,2)剪切模量GO(G)AEU(NM1,3)泊松比UO(u) AAI(NA,4)截面信息数组，其中： AAI(NA1,1)NA1截面面积A AAI(NA1,2)绕Z轴惯性矩IZ AAI(NA1,3)绕Y轴惯性矩IY AAI(NA1,4)绕X轴惯性矩IPP2(NP2,4)非节点荷载数组，其中：P2(I,1)作用节点号， P2(I,2)荷载类型码（见程序）P2(I,3)作用位置参数（见程序） P2(I,4)荷载大小JC(NC,2)约束信息数组，其中： JC(I,1)约束节点号 JC(I,2)约束代码取1，2，3，4，5，6代码 1:U=0, 2:V=0, 3:W=0, 4:SITA-X=0, 5:SITA-Y=0, 6:SITA-Z=0ANG(NE)单元角数组，即Y与Y所成的角度，用以求坐标变换矩阵AKEP(12,12)局部坐标系下单元刚度矩阵 AKE(12,12)整体坐标系下单元刚度矩阵 R(12,12)坐标变换矩阵 RT(12,12)R(12,12)的转置矩阵POP(12)局部坐标系下固端力向量 PO(12)整体坐标系下固端力向量A(NT)一维变带宽存储总刚度矩阵一维数组P(N)节点荷载向量 UPE(12)局部坐标系下节点位移IS(12)定位向量数组 LD(600)指示数组三、 源程序清单：C 3D-FRAME STRUCTURE ANALYAISPROGRAM MAIN /*主程序开始*/ CHARACTER*20 FNAME1,FNAME2 /*文件名变量说明*/ DIMENSION MEA(100,4),JZ(50,2),AAI(10,4),AEU(10,4),P2(50,4) /*数组说明*/ DIMENSION ANG(100) DIMENSION R(12,12),RT(12,12),RA(12,12),POP(12),PO(12) DIMENSION AKEP(12,12),AKE(12,12),P(600),X(100),Y(100),Z(100) DIMENSION LD(600),IS(12),A(50000),UPE(12),PPE(12),UE(12) WRITE(*,(A22)INPUT DATA FILE: READ(*,(A) FNAME1 /*输入原始数据文件名*/ WRITE(*,(A22)OUTPUT DATA FILE: READ(*,(A) FNAME2 /*输入结果数据文件名*/ WRITE(*,19) FNAME1,FNAME2 /*输出文件名以便较核*/ 19 FORMAT(1X,A,A) OPEN(1,FILE=FNAME1,STATUS=OLD) /*打开原始数据文件*/ OPEN(2,FILE=FNAME2,STATUS=NEW) /*打开结果数据文件*/ READ(1,*) NN,NE,NM,NA,NC,NP2 /*以下输入原始数据，无格式 NF=6 具体符号含义见前一页*/ ND=2 NFD=NF*ND N=NN*NF DO 50 I=1,NN READ(1,*) K,X(I),Y(I),Z(I),(P(NF*(I-1)+J),J=1,NF) 50 CONTINUE READ(1,*) (JZ(I,1),JZ(I,2),I=1,NC) DO 100 I=1,NE 100 READ(1,*) IE,(MEA(I,J),J=1,4) READ(1,*) (P2(I,J),J=1,4),I=1,NP2) DO 111 I=1,NE 111 READ(1,*) K,ANG(I) READ(1,*) (AAI(I,J),J=1,4),I=1,NA) READ(1,*) (AEU(I,J),J=1,3),I=1,NM) CLOSE(1) /*关闭原始数据文件*/ WRITE(2,460) NN,NE,NM,NA,NC,NP2 /*以下输出原始数据以便较核*/ 460 FORMAT(/3X,The input of the NN,NE,NM,NA,NC,NP2 is:/(5X,6I5) WRITE(2,465) (I,X(I),Y(I),Z(I),P(6*I-5),P(6*I-4), &P(6*I-3),P(6*I-2),P(6*I-1),P(6*I),I=1,NN) 465 FORMAT(/3X,The input of X(N),Y(N),Z(N),P(N) is:/(5X,I5,9F7.2) WRITE(2,470) (JZ(I,J),J=1,2),I=1,NC) 470 FORMAT(/3X,The input of the constrained information JZ is: & /(5X,I5,5X,I5) WRITE(2,475)(I,(MEA(I,J),J=1,4),I=1,NE) 475 FORMAT(/3X,The input of the element information MEA is: & /(5X,I5,5X,4I5) WRITE(2,480)(P2(I,J),J=1,4),I=1,NP2) 480 FORMAT(/3X,The input of the non-node load P2 is:/(5X,4F10.2) WRITE(2,485)(I,ANG(I),I=1,NE) 485 FORMAT(/3X,The element angle ANG is:/(5X,I5,F10.2) WRITE(2,490)(I,AAI(I,1),AAI(I,2),AAI(I,3),AAI(I,4),I=1,NA) 490 FORMAT(/3X,The input of the section information AAI is: & /(5X,I5,4(2X,E11.4) WRITE(2,495)(I,(AEU(I,J),J=1,3),I=1,NM) 495 FORMAT(/3X,The input of the materials information AEU is: & /(5X,I5,3(2X,E11.4)C STRUCTURE STIFFNESS STATEMENT CALL FLD(NN,NE,MEA,NF,ND,N,NT,LD) /*求LD数组*/DO 500 I=1,NT /*对A数组赋零*/ 500 A(I)=0.0 DO 600 IE=1,NE /*对单元进行循环以求结构刚度矩阵*/ CALL KEP(NN,IE,MEA,AEU,AAI,X,Y,Z,AKEP) /*局部坐标系下单元刚度矩阵*/ CALL CR(NN,IE,MEA,X,Y,Z,ANG,R) /*坐标转换矩阵*/ CALL TRAN(12,12,R,RT) CALL DOT(12,12,12,RT,AKEP,RA) CALL DOT(12,12,12,RA,R,AKE) /*整体坐标系下单元刚度矩阵*/ CALL FIS(IE,MEA,NF,ND,NFD,IS) / *定位向量数组IS*/ DO 560 I=1,NFD DO 560 J=1,NFD IF(IS(I)-IS(J) 560,520,520 /*整体单元刚度的下三角部分加到总刚中去*/ 520 NI=IS(I) IJ=LD(NI)-NI+IS(J) A(IJ)=A(IJ)+AKE(I,J) /*形成整体刚度矩阵 A*/ 560 CONTINUE 600 CONTINUEC EQUIVALENT NODE FORCE /*对单元循环求等效节点荷载*/ DO 700 IE=1,NE CALL FIXF(NN,IE,MEA,X,Y,Z,NP2,P2,POP) /*求局部坐标系下固端力向量*/ CALL CR(NN,IE,MEA,X,Y,Z,ANG,R) CALL TRAN(12,12,R,RT) CALL DOT(12,12,1,RT,POP,PO) /*整体坐标系下固端力向量*/ CALL FIS(IE,MEA,NF,ND,NFD,IS) DO 650 I=1,NFD NI=IS(I) P(NI)=P(NI)-PO(I) /*节点力叠加非节点力形成荷载向量*/ 650 CONTINUE 700 CONTINUE CALL FCC(NC,N,NT,NF,JZ,LD,A) /*处理位移约束条件*/ CALL DECOM(N,NT,A,P,LD) /*求解方程的节点位移*/ WRITE(2,850)(I,P(6*I-5),P(6*I-4),P(6*I-3),P(6*I-2),P(6*I-1), /*输出位移*/ &P(6*I),I=1,NN) 850 FORMAT(/3X,The nodal displacement is :/3X,Node NO., &2X,u,9X,v,9X,w,7X,Theta-x,5X,Theta-y,5X,Theta-z &/(3X,I4,6E11.4) WRITE(2,860) 860 FORMAT(/3X,The node inner force of element is:/1X,ELE NO., &I-J,4X,Fx,9X,Fy,9X,Fz,9X,Mx,9X,My,9X,Mz) DO 800 IE=1,NE /*对单元进行循环以求节点内力*/ CALL KEP(NN,IE,MEA,AEU,AAI,X,Y,Z,AKEP) /*局部坐标系下单元刚度矩阵*/ CALL CR(NN,IE,MEA,X,Y,Z,ANG,R) /*坐标转换矩阵*/ CALL FIS(IE,MEA,NF,ND,NFD,IS) /*定位向量*/ DO 750 I=1,NFD /*对单元自由度循环*/ NI=IS(I) UE(I)=P(NI) /*从整体位移中分离出单元位移*/ 750 CONTINUE CALL DOT(12,12,1,R,UE,UPE) /*局部坐标系下单元位移*/ CALL DOT(12,12,1,AKEP,UPE,PPE) /*局部坐标系下单元节点力*/ CALL FIXF(NN,IE,MEA,X,Y,Z,NP2,P2,POP) /*求固端力向量*/ DO 760 I=1,NFD PPE(I)=PPE(I)+POP(I) /*形成最后节点力向量*/ 760 CONTINUE WRITE(2,880) IE,(PPE(I),I=1,6) /*输出节点力*/ WRITE(2,900) (PPE(I),I=7,12) 800 CONTINUE 880 FORMAT(2X,I3,4X,I,6E11.4) 900 FORMAT(9X,J,6E11.4) CLOSE(2) /*关闭结果文件*/ END /*主程序结束*/C BEAM(3-D) ELEMENT LOCAL STIFNESS /*局部坐标系下单元刚度矩阵计算子程序*/ SUBROUTINE KEP(NN,IE,MEA,AEU,AAI,X,Y,Z,AKEP) DIMENSION X(NN),Y(NN),Z(NN),MEA(100,4),AEU(10,3),AAI(10,4) /*变量说明*/ DIMENSION AKEP(12,12) REAL IZ,IY,IP CALL ZERO(12,12,AKEP) /*对AKEP赋零*/ NM1=MEA(IE,4) /*IE单元的材料类型号*/ EO=AEU(NM1,1) /*弹性模量E*/ GO=AEU(NM1,2) /*剪切模量G*/ UO=AEU(NM1,3) /*泊松比*/ NA1=MEA(IE,3) /*单元IE截面类型号*/ AO=AAI(NA1,1) /*面积A*/ IZ=AAI(NA1,2) /*绕Z轴惯性矩IZ*/ IY=AAI(NA1,3) /*绕Y轴惯性矩IY*/ IP=AAI(NA1,4) /*绕X轴惯性矩IP*/ I=MEA(IE,1) /*节点 J=MEA(IE,2) 编号*/ SL=SQRT(X(J)-X(I)*2+(Y(J)-Y(I)*2+(Z(J)-Z(I)*2) /*IE单元梁长度*/ AKEP(1,1)=EO*AO/SL /*以下计算单元刚度矩阵元素*/ AKEP(7,1)=-AKEP(1,1) /*详细公式参阅温瑞槛所编 AKEP(7,7)=AKEP(1,1) 冶金工业出版社1998年出版*/ AKEP(2,2)=12*EO*IZ/(SL*3.0) AKEP(8,2)=-AKEP(2,2) AKEP(8,8)=AKEP(2,2) AKEP(6,2)=6*EO*IZ/(SL*2) AKEP(12,2)=AKEP(6,2) AKEP(8,6)=-AKEP(6,2) AKEP(12,8)=-AKEP(6,2) AKEP(3,3)=12*EO*IY/(SL*3) AKEP(9,3)=-AKEP(3,3) AKEP(9,9)=AKEP(3,3) AKEP(5,3)=-6*EO*IY/(SL*2) AKEP(11,3)=AKEP(5,3) AKEP(9,5)=-AKEP(5,3) AKEP(11,9)=-AKEP(5,3) AKEP(4,4)=GO*IP/SL AKEP(10,4)=-AKEP(4,4) AKEP(10,10)=AKEP(4,4) AKEP(5,5)=4*EO*IY/SL AKEP(11,5)=AKEP(5,5)/2.0 AKEP(11,11)=AKEP(5,5) AKEP(6,6)=4*EO*IZ/SL AKEP(12,6)=AKEP(6,6)/2.0 AKEP(12,12)=AKEP(6,6) DO 10 II=1,11 /*AKEP为对称矩阵，循环求上三角*/ DO 10 JJ=II+1,12 AKEP(II,JJ)=AKEP(JJ,II) 10 CONTINUE RETURN ENDC BEAM(3-D) ELEMENT COORDINATE TRANSFORMATION /*坐标转换矩阵计算子程序*/ SUBROUTINE CR(NN,IE,MEA,X,Y,Z,ANG,R) DIMENSION X(NN),Y(NN),Z(NN),MEA(100,4),ANG(100),R(12,12) /*变量说明*/ REAL LMD CALL ZERO(12,12,R) /*初始赋零*/ I=MEA(IE,1) /*节点编号*/ J=MEA(IE,2) XL=X(J)-X(I) YL=Y(J)-Y(I) ZL=Z(J)-Z(I) SL=SQRT(XL*XL+YL*YL+ZL*ZL) /*IE单元梁长度*/ COX=XL/SL /*三个轴的方向角*/ COY=YL/SL COZ=ZL/SL LMD=SQRT( COX*COX+COY*COY) ANG1=ANG(IE) CO=COS(ANG1) SI=SIN(ANG1) IF(LMD.EQ.0.0) THEN /*竖直杆单元*/ R(1,3)=1.0 R(2,1)=-SI R(2,2)=CO R(3,1)=-CO R(3,2)=-SI ELSE /*一般杆单元*/ R(1,1)=COX /*详细公式参阅温瑞槛所编 R(1,2)=COY 冶金工业出版社1998年出版*/ R(1,3)=COZ R(2,1)=-COY/LMD R(2,2)=COX/LMD R(2,3)=0.0 R(3,1)=-COX*COZ/LMD R(3,2)=-COY*COZ/LMD R(3,3)=LMD END IF DO 70 II=1,3 /*循环把R冲满*/ DO 60 JJ=1,3 R(II+3,JJ+3)=R(II,JJ) R(II+6,JJ+6)=R(II,JJ) R(II+9,JJ+9)=R(II,JJ) 60 CONTINUE 70 CONTINUE RETURN ENDC THE FOLLOWINF PARAGRAPH ZERO A MATRIX /*矩阵赋零子程序*/ SUBROUTINE ZERO(M,N,A) DIMENSION A(M,N) DO 20 II=1,M DO 20 JJ=1,N A(II,JJ)=0.0 20 CONTINUE RETURN ENDC MATRIX TRANSFORMATION /*求矩阵转置子程序*/ SUBROUTINE TRAN(M,N,A,B) DIMENSION A(M,N),B(N,M) DO 30 II=1,M DO 30 JJ=1,N B(JJ,II)=A(II,JJ) 30 CONTINUE RETURN ENDC MATRIX MULTIPLY /*矩阵相乘子程序*/ SUBROUTINE DOT(M,N,L,A,B,C) DIMENSION A(M,N),B(N,L),C(M,L) DO 50 I=1,M DO 50 J=1,L C(I,J)=0.0 DO 40 K=1,N C(I,J)=C(I,J)+A(I,K)*B(K,J) 40 CONTINUE 50 CONTINUE RETURN ENDC FIXED END FORCE /*固端反力求解子程序*/ SUBROUTINE FIXF(NN,IE,MEA,X,Y,Z,NP2,P2,POP) DIMENSION X(NN),Y(NN),Z(NN),P2(50,4),MEA(100,4),POP(12) /*变量说明*/ I=MEA(IE,1) /*求节点编号*/ J=MEA(IE,2) SL=SQRT(X(J)-X(I)*2+(Y(J)-Y(I)*2+(Z(J)-Z(I)*2) /*单元长度*/ CALL ZERO(12,1,POP) DO 110 II=1,NP2 /*对每一个非节点力进行循环*/ IF(IE.EQ.P2(II,1) THEN /*看II非节点力是否作用在所求单元上*/ IF(P2(II,2).EQ.1.0).OR.(P2(II,2).EQ.2.0) THEN /*II非节点力为集中力*/ A=P2(II,3) /*作用位置*/ B=SL-A PP=P2(II,4) /*力的大小*/ PQI=-PP*(B*2)*(SL+2*A)/(SL*3) PQJ=-PP*(A*2)*(SL+2*B)/(SL*3) PMI=-PP*A*(B*2)/(SL*2) PMJ=PP*B*(A*2)/(SL*2) IF(P2(II,2).EQ.1.0) THEN /*1：集中力作用在XOY平面内*/ POP(2)=POP(2)+PQI POP(6)=POP(6)+PMI POP(8)=POP(8)+PQJ POP(12)=POP(12)+PMJ ELSE /*2：集中力作用在XOZ平面内*/ POP(3)=POP(3)+PQI POP(5)=POP(5)-PMI POP(9)=POP(9)+PQJ POP(11)=POP(11)-PMJ END IF ELSE IF(P2(II,2).EQ.3.0).OR.(P2(II,2).EQ.4.0) THEN /*荷载为均布力*/ A=P2(II,3) /*从X=0开始的作用长度*/ PQ=P2(II,4) /*集度（大小）*/ PQI=-PQ*A*(SL*3-SL*(A*2)+0.5*(A*3)/(SL*3) PQJ=-PQ*(A*3)*(1-A/(2*SL)/(SL*2) PMI=-PQ*(A*2)*(6-8*A/SL+3*A*A/(SL*SL)/12.0 PMJ=PQ*(A*3)*(4-3*A/SL)/(12.0*SL) IF(P2(II,2).EQ.3.0) THEN /*3：均布力作用在XOY平面内*/ POP(2)=POP(2)+PQI POP(6)=POP(6)+PMI POP(8)=POP(8)+PQJ POP(12)=POP(12)+PMJ ELSE /*4：均布力作用在XOZ平面内*/ POP(3)=POP(3)+PQI POP(5)=POP(5)-PMI POP(9)=POP(9)+PQJ POP(11)=POP(11)-PMJ END IF ELSE /*荷载为集中力偶矩（5，6）*/ A=P2(II,3) /*作用位置*/ B=SL-A PM=P2(II,4) /*大小*/ PQI=6*A*B*PM/(SL*3) PQJ=-PQI PMI=B*(3*A-SL)*PM/(SL*SL) PMJ=A*(3*B-SL)*PM/(SL*SL) IF(P2(II,2).EQ.5.0) THEN /*5：集中力偶矩作用在XOY平面内*/ POP(2)=POP(2)+PQI POP(6)=POP(6)+PMI POP(8)=POP(8)+PQJ POP(12)=POP(12)+PMJ ELSE /*6：集中力偶矩作用在XOZ平面内* POP(3)=POP(3)-PQI POP(5)=POP(5)+PMI POP(9)=POP(9)-PQJ POP(11)=POP(11)+PMJ END IF END IF END IF 110 CONTINUE RETURN ENDC IS DIMENSION /*单元位移与结构位移对应的IS数组计算子程序*/ SUBROUTINE FIS(IE,MEA,NF,ND,NFD,IS) DIMENSION MEA(100,4),IS(NFD) DO 150 ID=1,ND /*对单元节点循环*/ DO 120 JF=1,NF /*对节点自由度进行循环*/ I1=(ID-1)*NF+JF IS(I1)=(MEA(IE,ID)-1)*NF+JF /*计算IS数组的值*/ 120 CONTINUE 150 CONTINUE RETURN ENDC LD DIMENSION /*LD数组计算子程序*/ SUBROUTINE FLD(NN,NE,MEA,NF,ND,N,NT,LD) DIMENSION MEA(100,4),LD(N) LD(1)=1 DO 300 K=1,NN /*对节点进行循环查找NMIN节点*/ NMIN=5000 /*对NMIN赋初值，宜大*/ DO 200 I=1,NE /*对单元循环，判断是否有K点*/ DO 150 J=1,ND IF(MEA(I,J).NE.K) GOTO 150 DO 100 L=1,ND IF(MEA(I,L).LT.NMIN) NMIN=MEA(I,L) /*找与K有关的最小节点*/ 100 CONTINUE 150 CONTINUE 200 CONTINUE DO 250 I=1,NF J=(K-1)*NF+I IF(J.NE.1) LD(J)=LD(J-1)+(K-NMIN)*NF+I /*计算K点对应各行的LD值*/ 250 CONTINUE 300 CONTINUE NT=LD(N) /*A的容量*/ RETURN ENDC CONSTRITED CONDITION /*约束处理子程序*/ SUBROUTINE FCC(NC,N,NT,NF,JZ,LD,A) DIMENSION JZ(50,2),LD(N),A(NT) DO 350 K=1,NC /*对没个约束循环*/ I=JZ(K,1) /*约束的节点号*/ J=JZ(K,2) /*约束的自由度号*/ NI=(I-1)*NF+J NJ=LD(NI) A(NJ)=1.0E25 /*主对角元素置大值*/ 350 CONTINUE RETURN ENDC MATIX DECOMPOSITION FUNCTION /*矩阵分解法解线性方程组子程序 SUBROUTINE DECOM(N,NT,A,B,LD) 详细注解请参考李章政所编写 DIMENSION A(NT),B(N),LD(N) 结构分析程序设计2003年2月第二版*/ DO 20 I=1,N LDI=LD(I) IF(I.EQ.1) GOTO 10 IO=LDI-I IM1=I-1 MI=1-IO+LD(IM1) IF(MI.GT.IM1) GOTO 10 DO 8 J=MI,IM1 JO=LD(J)-J IJ=IO+J IF(J.EQ.1) GOTO 6 JM1=J-1 MJ=1-JO+LD(JM1) MIJ=MAX0(MI,MJ) IF(MIJ.GT.JM1) GOTO 6 S=0.0 DO 2 K=MIJ,JM1