有限元计算结构力学fortran程序.docx

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

2、! E弹性模量； ! JC单元定位向量数组； ! 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、正负; ! 4三角形荷载，c为荷载距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

5、(I,2)-类型；PF(I,3)-误差大小(含正负),正负取决于消除 ! 误差时端点的运动方向，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中的

6、地址，JI=KD(J)-J+I； ! 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),& K

7、D(150),TK(1000),P(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 &

8、OF STRUCTURE/5X,NE=,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

9、(TK,NN,KD,N) !对总刚进行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)

10、!解方程组并输出结点位移，存放在数组P(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,

11、3X,V2=,D12.4,3X,M2=,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

12、 READ(1,70) TL READ(1,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),

13、I=1,NE) WRITE(2,30)(I,(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 SPECIAL JOINTS/6X,JOINT,4X,u,4x,v,4x,ceta/(6X,4I5)70