计算结构力学程序.doc

《计算结构力学程序.doc》由会员分享，可在线阅读，更多相关《计算结构力学程序.doc（18页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、计算结构力学编程大作业时间： 2007年6月 !* ! 关于程序的说明 !* !一、功能： ! 1、可计算包括节点力，一般非节点力，支座沉降、温度荷载作用、制造误差的平 ! 面桁架、梁、刚架及其组合结构的节点位移与杆端力； ! 2、可同时计算多种工况下的节点位移与杆端力。 !* !* ! ! 二、变量说明： ! NE单元数； ! N结构中自由度数； ! NJ节点数； ! NS特殊节点数，包括支座节点、主从节点（1节点不做主节点）、连接桁架的铰节点（没有转角）； ! NAI结构的单元截面类型数； ! MT单元截面类型号； ! NL荷载工况数； ! H截面高度； ! E弹性模量； ! JC单元定位

2、向量数组； ! X(NJ),Y(NJ)节点的X，Y坐标值； ! JE(NE,2)单元两端节点码数组； ! AI(NAI,2)按单元类型顺序存放A与I，AI（I，1）第I类单元的截面积，AI（I，2）第I类单元的 ! 惯性矩； ! MT(NE)单元所属单元类型号； ! JS(NS,4)特殊节点信息，JS(I,1)结点码；JS(I,2)，JS(I,3)，JS(I,4)U，V，CETA约束信息， ! 有约束为1，没有约束为0；从节点某位移同主节点时位移时，该位移约束信息填主节点码； ! ! PJ(NP,3)节点荷载信息数组；PJ（I，1）节点力所在节点号；PJ（I，2）节点力作用坐标方向： ! 坐标

3、方向U，V，M分别为1，2，3； PJ（I，3）节点力的大小(含正负号)；U，V方向集中力时， ! 与坐标轴正向同向为正，M按右手法则为正；本程序推导过程取y轴向下为正。 ! ! PF(NF,4)非节点荷载数组，并给出以下类型说明： ! 前6类型数据输法(梯形等可以用叠加法计算)： ! PF(I,1)-单元码；PF(I,2)-类型；PF(I,3)-荷载大小；PF(I,4)-c值； ! 1垂直于单元的均布力,大小为q，以坐标轴正向为正,c为荷载末端距i节点距离; ! 2非节点集中力P，c为荷载距i节点距离; ! 3非节点集中力距M，c为荷载距i节点距离，右手法则判正负; ! 4三角形荷载，c为荷

4、载距i节点距离，i端为0，距离i端c时力为q; ! j端为0的三角形，可按叠加法处理。 ! 5沿杆轴向均布力,大小为q,c为荷载末端距i节点距离; ! 6沿杆轴向集中力,大小为q,c为荷载末端距i节点距离; ! ! 从第7到第9类型（支座沉降）数据输法：PF(I,1)-单元码；PF(I,2)-类型；PF(I,3)-位移大小（含正负），坐标轴正向位为正，转角按右手法则；PF(I,4)-沉降所在的单元位移分量，i端为1-3，j端为4-6； ! ! 7沿轴向支座沉降； ! 8垂直于轴向支座沉降； ! 9支座转动； ! 10制造误差,PF(I,1)制造误差所在单元，PF(I,2)-类型；PF(I,3)

5、-误差大小(含正负),正负取决于消除 ! 误差时端点的运动方向，PF(I,4)误差所在坐标号； ! 11温度荷载，PF(I,1)荷载所在单元，数据形式为：ElementNo.1 ,如2单元上有温度荷载，则PF(I,1)=2.1; ! PF(I,2)温度变化值t1, PF(I,3)温度变化值t2,PF(I,4)材料线膨胀系数; ! ! TK(NN)采用一维存储结构刚度矩阵，上半带元素（每列第一个非零元素到对角元）； ! KD主元地址数组，表示结构刚度矩阵的主元在TK中的序号，KD中最后一个数是TK中元素的总个数； ! JI结构刚度矩阵上半带的非对角元素在TK中的地址，JI=KD(J)-J+I；

6、! JN(NJ,3)结点位移分量编号数组，用于存放结点三个位移的位移分量号码， ! JN(I,1)，JN(I,2)，JN(I,2)-分别为结点I的U，V，CETA分量的位移分量（坐标）号码； ! ! P(N)节点荷载列阵；在回代求位移时存放位移量； ! F(N)求得的杆端力列阵； ! FO(6)等效节点荷载列阵； ! ! !* !* 平面结构分析源程序内容 * !* PROGRAM PFF DIMENSION X(50),Y(50),JE(50,2),MT(30),AI(10,2),JS(20,4),PJ(50,3),PF(50,4),JN(50,3),& KD(150),TK(1000),P

7、(150),F(6),H(50) DOUBLE PRECISION TK,P,F CHARACTER *200 TL OPEN(1,FILE=INDAT.DAT,STATUS=OLD) OPEN(2,FILE=OUTDAT.DAT,STATUS=NEW) READ(1,70) TL READ(1,70) TL READ(1,*)NE,NJ,NS,NAI,NL,E WRITE(2,10)NE,NJ,NS,NAI,NL,E10 FORMAT(5X,PLANE FRAME STRUCTURE ANALYSIS/5X,*/2X,CONTROL PARAMETERS &OF STRUCTURE/5X,N

8、E=,I2,8X,NJ=,I2,8X,NS=,I2,8X,NAI=,I2,/5X,NL=,I2,8X,E=,E12.4) CALL INPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H) !读入数据文件 CALL DJN(NJ,NS,JS,JN,N) !计算结构自由度数N，形成结点位移分量数组JN CALL ADE(NE,NJ,N,JE,JN,KD,NN)!形成主元地址数组KD(N) CALL SSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK) !形成总刚，一维存储数组TK(NN) CALL UTDU3(TK,NN,KD,N) !对总刚进

9、行UTDU分解，以用于解方程组 DO 20 LC=1,NL !对工况循环 READ(1,70)TL READ(1,70)TL READ(1,70)TL READ(1,*)NP,NF !读入工况信息 WRITE(2,30)LC,NP,NF30 FORMAT(/2X,LOAD DATA/10X,LOAD CASE=,I3/10X,NP=,I3,8X,NF=,I3) CALL NLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H) !形成总荷载列阵P(N) CALL BACK3(TK,NN,P,N,KD,JN,NJ)!解方程组并输出结点位移，存放在数组P

10、(N)中 WRITE(2,40)40 FORMAT(/4X,MEMBER-END FORCES OF ELEMENTS/4X,ELEMENT,13X,N,17X,V,17X,M) DO 60 M=1,NE CALL MQN(M,NE,NJ,NAI,N,NF,E,X,Y,JE,MT,AI,JN,PF,P,F,H)!计算单元杆端力，存放在数组F(6)中 WRITE(2,50)M,(F(I),I=1,6) !输出杆端力50 FORMAT(/1X,I10,3X,N1=,D12.4,3X,V1=,D12.4,3X,M1=,D12.4/14X,N2=,&D12.4,3X,V2=,D12.4,3X,M2=,

11、D12.4)60 CONTINUE20 CONTINUE70 FORMAT(A) CLOSE(1) CLOSE(2) END SUBROUTINE INPUT(NE,NJ,NS,NAI,X,Y,JE,MT,AI,JS,H)!读入数据文件 DIMENSION X(NJ),Y(NJ),JE(NE,2),MT(NE),AI(NAI,2),JS(NS,4),H(NE) INTEGER NO READ(1,70) TL READ(1,70) TL READ(1,70) TL READ(1,*)(NO,X(I),Y(I),I=1,NJ) READ(1,70) TL READ(1,70) TL READ(1

12、,70) TL READ(1,*)(NO,JE(I,1),JE(I,2),MT(I),H(I),I=1,NE) READ(1,70) TL READ(1,70) TL READ(1,70) TL READ(1,*)(NO,(AI(I,J),J=1,2),I=1,NAI) READ(1,70) TL READ(1,70) TL READ(1,70) TL READ(1,*)(JS(I,J),J=1,4),I=1,NS) WRITE(2,10)(I,X(I),Y(I),I=1,NJ) WRITE(2,20)(I,JE(I,1),JE(I,2),MT(I),I=1,NE) WRITE(2,30)(I

13、,(AI(I,J),J=1,2),I=1,NAI) WRITE(2,40)(JS(I,J),J=1,4),I=1,NS)10 FORMAT(/2X,COORDINATES OF JOINTS/6X,JOINT,11X,X,11X,Y/(6X,I4,5X,2F12.4)20 FORMAT(/2X,INFORMATION OF ELEMENTS/6X,ELEMENT,& 4X,JOINT-I,4X,JOINT-J,5X,TYPE/(2X,4I10)30 FORMAT(/7X,TYPE,10X,A,12X,I/(8X,I2,5X,2F12.6)40 FORMAT(/2X,INFORMATION OF

14、 SPECIAL JOINTS/6X,JOINT,4X,u,4x,v,4x,ceta/(6X,4I5)70 FORMAT(A) END! Determine number of joint displacements. SUBROUTINE DJN(NJ,NS,JS,JN,N) !计算结构自由度数N，形成数组JN。JN为结点位移分量编号数组，用于存放结点三个位移的位移分量号码 DIMENSION JS(NS,4),JN(NJ,3) DO 10 I=1,NJ DO 10 J=1,310 JN(I,J)=0 !对JN置0 DO 20 I=1,NS L=JS(I,1) !取特殊节点的节点码 DO 2

15、0 J=1,320 JN(L,J)=JS(I,J+1) !将JS中位移约束信息JS(I,2),JS(I,3),JS(I,4)附给JN(L,1),JN(L,2),JN(L,3) N=0 !自由度置0 ID=0 DO 30 I=1,NJ DO 30 J=1,3 IF(JN(I,J)-1)40,50,60 !由条件的小于0，等于0，大于0来执行40，50，60语句40 N=N+1 !小于0时，JN(I,J)=0，有自由度，自由度累加 JN(I,J)=N GOTO 3050 JN(I,J)=0 !等于0，JN(I,J)=1，没有自由度 GOTO 3060 ID=1 !大于0，结点I是从结点,自由度不增

16、加，RETURN30 CONTINUE IF(ID.EQ.0)GOTO 80 !判断是否有主从结点 DO 70 I=1,NS L=JS(I,1) DO 70 J=1,3 K=JS(I,J+1) IF(K.LE.1)GOTO 70 !若K1，则L是从结点，K是主结点编号 JN(L,J)=JN(K,J) !从结点L的第J个位移分量的编号应取主结点K的第J个位移分量的编号70 CONTINUE80 RETURN END! Determine address of diagonal elements of total stiffness matrix. SUBROUTINE ADE(NE,NJ,N,J

17、E,JN,KD,NN) !形成主元地址数组KD(N) DIMENSION JE(NE,2),JN(NJ,3),KD(N),JC(6) DO 10 I=1,N10 KD(I)=0 !主元地址数组KD(N)先置0 DO 20 M=1,NE !对每根杆件循环 CALL EJC(M,NE,NJ,JE,JN,JC) !形成单元定位向量JC MIN=N !单元定位向量中的最小非零分量MIN，赋初值为N DO 30 I=1,6 !对定位向量的六个分量循环 J=JC(I) !取定位向量的分量 IF(J.EQ.0) GOTO 30 !排除零分量 IF(J.LT.MIN) MIN=J !求定位向量的最小非零分量3

18、0 CONTINUE DO 20 I=1,6 !对定位向量的六个分量循环 J=JC(I) !取定位向量的分量 IF(J.EQ.0) GOTO 20 !排除零分量 NW=J-MIN+1 !计算半带宽 IF(NW.GT.KD(J) KD(J)=NW !第J列的真正半带宽20 CONTINUE NN=1!NN置初值为1，NN为存放总刚度矩阵中非零元素的一维数组TK(NN)的元素个数 DO 40 J=2,N !对总刚度矩阵K列循环 NN=NN+KD(J) KD(J)=NN !累加各列半带宽，形成KD数组40 CONTINUE END! Assemble structure stiffness matr

19、ix stored as a vector. SUBROUTINE SSM(NE,NJ,NAI,E,N,NN,X,Y,JE,MT,AI,JN,KD,TK) !形成存储总刚的一维数组TK(NN) DIMENSION X(NJ),Y(NJ),JE(NE,2),JN(NJ,3),MT(NE),AI(NAI,2),KD(N),TK(NN),KE(6,6),JC(6) DOUBLE PRECISION TK,KE,BL,SI,CO DO 10 I=1,NN10 TK(I)=0.0 !数组TK(NN)置0 DO 20 M=1,NE !对每根杆件循环 CALL SCL(M,NE,NJ,X,Y,JE,BL,S

20、I,CO) !求单元常数 CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成global坐标下的单元刚度矩阵KE(6,6) CALL EJC(M,NE,NJ,JE,JN,JC) !形成单元定位向量JC DO 30 L=1,6 !对单元刚度矩阵的行循环 I=JC(L) !取单刚的第L行在总刚中的行码 IF(I.EQ.0) GOTO 30 DO 40 K=1,6 !对单元刚度矩阵的列循环 J=JC(K) !取单刚的第K列在总刚中的列码 IF(J.LT.I) GOTO 40 !排除总刚的下三角元素 JI=KD(J)-J+I !计算总刚第I行第J列元素在数组TK中的地

21、址 TK(JI)=TK(JI)+KE(L,K) !叠加单刚元素到总刚中40 CONTINUE30 CONTINUE20 CONTINUE ! WRITE(2,50) TK!50 FORMAT(TK,3X,2D15.4) END SUBROUTINE SCL (M,NE,NJ,X,Y,JE,BL,SI,CO) !计算单元常数 DIMENSION X(NJ),Y(NJ),JE(NE,2) DOUBLE PRECISION BL,SI,CO I=JE(M,1) !取M单元起始结点号 J=JE(M,2) !取M单元终止结点号 DX=X(J)-X(I) DY=Y(J)-Y(I) BL=SQRT(DX*D

22、X+DY*DY)!计算单元杆长 SI=DY/BL !计算单元转角SIN值 CO=DX/BL !计算单元转角COS值 END! Decompose total stiffness matrix using Crout method. SUBROUTINE UTDU3(A,NN,ID,N)!对总刚进行UTDU分解，用于解方程组，对应实参(TK,NN,KD,N) DIMENSION A(NN),ID(N) DOUBLE PRECISION A DO 30 J=2,N !对列循环 JJ=ID(J) !取主对角元素在一维数组A(NN)中的地址码,A(NN)对应TK(NN) J1=J-1 MJ=J-JJ+

23、ID(J1)+1 !第J列第一个非零元素的行号 IF(MJ.GT.J1) GOTO 30 !排除下三角元素 DO 20 I=MJ,J !对行循环 II=ID(I) I1=I-1 MIJ=MJ IF(I.EQ.J) GOTO 5 MI=I-II+ID(I1)+1 IF(MIJ.LT.MI) MIJ=MI5 S=0.0 IF(MIJ.GT.I1) GOTO 15 DO 10 K=MIJ,I1 KI=II-I+K KK=ID(K) KJ=JJ-J+K10 S=S+A(KI)*A(KK)*A(KJ)15 IF(I.EQ.J) THEN A(JJ)=A(JJ)-S !求对角阵D的第J个元素，存放在一维数

24、组A(NN)中 ELSE JI=JJ-J+I A(JI)=(A(JI)-S)/A(II) !求三角阵U的第I行第J列元素，存放在一维数组A(NN)中 END IF20 CONTINUE30 CONTINUE END! Form total joint load vector. SUBROUTINE NLV(NE,NJ,NAI,E,N,NP,NF,X,Y,JE,JN,PJ,PF,MT,AI,P,H) !形成总荷载列阵P(N) DIMENSION X(NJ),Y(NJ),JE(NE,2),JN(NJ,3),PJ(NP,3),PF(NF,4),MT(NE),AI(NAI,2),& P(N),FO(6

25、),FE(6),JC(6),H(NE) DOUBLE PRECISION P,FO,FE,BL,SI,CO INTEGER NO DO 10 I=1,N P(I)=0.0 !总荷载列阵P(N)置010 CONTINUE READ(1,70)TL READ(1,70)TL READ(1,70)TL IF(NP.GT.0) THEN READ(1,*)(NO,(PJ(I,J),J=1,3),I=1,NP) !读入结点荷载信息 WRITE(2,20)(PJ(I,J),J=1,3),I=1,NP)20 FORMAT(/2X,JOINT LOAD/6X,JOINT,8X,XYM,12X,LOAD/(6X

26、,F5.0,6X,F5.0,6X,F12.4) DO 30 I=1,NP J=INT(PJ(I,1) !取结点荷载所在的结点号 K=INT(PJ(I,2) !取结点荷载在结点上的方向 L=JN(J,K) !取结点荷载的位移编号 IF(L.NE.0) P(L)=PJ(I,3) !把结点荷载叠加到总荷载列阵P(N)中30 CONTINUE END IF READ(1,70)TL READ(1,70)TL READ(1,70)TL IF(NF.GT.0) THEN READ(1,*)(NO,(PF(I,J),J=1,4),I=1,NF) !读入非结点荷载信息 WRITE(2,40)(PF(I,J),

27、J=1,4),I=1,NF)40FORMAT(/2X,NON-JOINT LOAD/6X,ELEMENT,8X,TYPE,8X,LOAD,&12X,C/(6X,F6.0,6X,F6.0,4X,2F12.4) DO 50 I=1,NF !对非结点荷载循环 M=INT(PF(I,1) !取非结点荷载所在的单元码 CALL SCL(M,NE,NJ,X,Y,JE,BL,SI,CO) !计算单元常数 CALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H)!将非结点荷载转化为等效结点荷载，存放在FO(6)中 CALL EJC(M,NE,NJ,JE,JN,JC)!形成单元定位数

28、组 FE(1)=-FO(1)*CO+FO(2)*SI!将等效结点荷载转化到global坐标下,y轴向下的坐标系。 FE(2)=-FO(1)*SI-FO(2)*CO!将等效结点荷载转化到global坐标下 FE(3)=-FO(3)!将等效结点荷载转化到global坐标下 FE(4)=-FO(4)*CO+FO(5)*SI!将等效结点荷载转化到global坐标下 FE(5)=-FO(4)*SI-FO(5)*CO!将等效结点荷载转化到global坐标下 FE(6)=-FO(6)!将等效结点荷载转化到global坐标下 DO 60 J=1,6 L=JC(J) IF(L.NE.0) P(L)=P(L)+F

29、E(J)!将等效结点荷载叠加到总荷载列阵P(N)中60 CONTINUE50 CONTINUE END IF70 FORMAT(A) END! Solve equations and print joint displacements. SUBROUTINE BACK3(A,NN,B,N,ID,JN,NJ) !解方程组并输出结点位移，存放在数组B(N)中。这里B(N)数组对应主程序中的P(N)数组 DIMENSION A(NN),B(N),ID(N),JN(NJ,3),D(3) DOUBLE PRECISION A,D,B DO 20 J=2,N!对列循环 JJ=ID(J) J1=J-1 MJ

30、=J-JJ+ID(J1)+1 !第J列第一个非零元素的行号 IF(MJ.GT.J1) GOTO 20 !若MJJ-1，则z的第J个值即为B的第J个值 DO 10 I=MJ,J1 JI=JJ-J+I10 B(J)=B(J)-A(JI)*B(I) !求z的第J个值，存放在B(N)数组。这里B(N)数组对应主程序中的P(N)数组20 CONTINUE DO 30 I=1,N II=ID(I) B(I)=B(I)/A(II) !求y的第J个值，存放在B(N)数组30 CONTINUE NW=1 !半带宽置初值为1 DO 35 I=2,N IW=ID(I)-ID(I-1) IF(IW.GT.NW) NW

31、=IW !计算最大列半带宽35 CONTINUE N1=N-1 DO 50 I=N1,1,-1 !对行循环 MIN=I+NW-1 IF(N.LT.MIN) MIN=N DO 40 J=I+1,MIN !对列循环 JJ=ID(J) MJ=J-JJ+ID(J-1)+1 !第J列第一个非零元素的行号 IF(I.LT.MJ) GOTO 40 !判断是否为带内元素 JI=JJ-J+I B(I)=B(I)-A(JI)*B(J) !求x的第I个值，存放在B(N)数组。这里x的值即为结点位移40 CONTINUE50 CONTINUE WRITE(2,60)60 FORMAT(/2X,JOINT DISPLA

32、CEMENTS/5X,JOINT,12X,U,14X,V,11X,ceta) DO 80 I=1,NJ DO 70 J=1,3 D(J)=0.0 L=JN(I,J) !取位移号 IF(L.NE.0) D(J)=B(L) !给位移赋值70 CONTINUE WRITE(2,75) I,D(1),D(2),D(3) !输出第I号结点的位移值75 FORMAT(2X,I6,4X,3D15.6)80 CONTINUE END SUBROUTINE MQN(M,NE,NJ,NAI,N,NF,E,X,Y,JE,MT,AI,JN,PF,P,F,H)!计算单元杆端力，存放在数组F(6)中 DIMENSION

33、X(NJ),Y(NJ),JE(NE,2),AI(NAI,2),P(N),JN(NJ,3),PF(NF,4),JC(6),FO(6),D(6),FD(6),&F(6),KE(6,6),MT(NE),H(NE) DOUBLE PRECISION P,FO,D,FD,F,KE,BL,SI,CO CALL SCL(M,NE,NJ,X,Y,JE,BL,SI,CO) !求单元常数 CALL ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE)!形成global坐标下的单元刚度矩阵KE(6,6) CALL EJC(M,NE,NJ,JE,JN,JC)!形成单元定位向量JC DO 10 I=1,6

34、 L=JC(I) !取位移编号 D(I)=0.0 IF(L.NE.0) D(I)=P(L) !取global坐标下的杆端位移10 CONTINUE DO 20 I=1,6 FD(I)=0.0 DO 30 J=1,6 FD(I)=FD(I)+KE(I,J)*D(J) !计算杆端力30 CONTINUE20 CONTINUE F(1)=FD(1)*CO+FD(2)*SI !将杆端力转化到local坐标下 F(2)=-FD(1)*SI+FD(2)*CO !将杆端力转化到local坐标下 F(3)=FD(3) !将杆端力转化到local坐标下 F(4)=FD(4)*CO+FD(5)*SI !将杆端力转

35、化到local坐标下 F(5)=-FD(4)*SI+FD(5)*CO !将杆端力转化到local坐标下 F(6)=FD(6) !将杆端力转化到local坐标下 IF(NF.GT.0)THEN DO 40 I=1,NF L=INT(PF(I,1) IF(M.EQ.L)THEN CALL EFX(I,M,NE,NAI,E,NF,PF,MT,AI,BL,FO,H) DO 50 J=1,6 F(J)=F(J)+FO(J) !如果有非结点荷载，杆端力中要加入等效结点荷载50 CONTINUE ENDIF40 CONTINUE ENDIF END SUBROUTINE EJC(M,NE,NJ,JE,JN,

36、JC) !形成单元定位向量JC，JC存放的是单元两端点的位移号 DIMENSION JE(NE,2),JN(NJ,3),JC(6) I=JE(M,1) !取M单元起始结点号 J=JE(M,2) !取M单元终止结点号 DO 10 K=1,3 JC(K)=JN(I,K) !取单元位移号赋值给定位向量JC JC(K+3)=JN(J,K) !取单元位移号赋值给定位向量JC10 CONTINUE END SUBROUTINE ESM(M,NE,NAI,E,MT,AI,BL,SI,CO,KE) !形成global坐标下的单元刚度矩阵KE(6,6) DIMENSION MT(NE),AI(NAI,2),KE

37、(6,6) DOUBLE PRECISION BL,SI,CO,KE I=MT(M) !取M单元截面类型号，MT单元截面类型号 C1=E*AI(I,1)/BL C2=2.0*E*AI(I,2)/BL C3=3.0*C2/BL C4=2.0*C3/BL S1=C1*CO*CO+C4*SI*SI S2=(C1-C4)*SI*CO S3=C3*SI S4=C1*SI*SI+C4*CO*CO S5=C3*CO S6=C2 KE(1,1)=S1 KE(1,2)=S2 KE(1,3)=-S3 KE(1,4)=-S1 KE(1,5)=-S2 KE(1,6)=-S3 KE(2,2)=S4 KE(2,3)=S5 KE(2,4)=-S2 KE(2,5)=-S4 KE(2,6)=S5 KE(3,3)=2.0*S6 KE(3,4)=S3 KE(3,5)=-S5 KE(3,6)=S6 KE(4,4)=S1 KE(4,5)=S2 KE(4,6)=S3 KE(5,5)=S4 KE(5,6)=-S5 KE(6,6)=2.