结构分析程序设计基础--论文.docx

结构分析程序设计基础姓 名： 学 号： 专 业： 授课教师： 指导教师： 2016年5月20日23一、结合自身所学专业，是阐述利用有限单元法进行结构分析时所常用离散模型及其特点。应用有限单元法求解土木工程问题时常用的离散模型有空间梁单元、板壳元、三维实体元及梁格单元。空间梁单元。空间梁单元是一种常用的单元,除了桥梁上部结构外,框架系统和动力联动装置都可以采用梁单元模拟。对于一个采用梁单元模拟的构件,一般它在一个方向上的尺寸最少是另两个方向上的10倍。梁单元除了承担拉力和压力，还可以承担剪力和弯矩。三维梁单元每个节点有6个自由度,即沿轴的位移和绕这3个轴的转角。根据结构受载后截面是否保持平截面，可分为自由扭转和翘曲扭转。后者考虑了受载后横截面双力矩阵和翘曲扭矩两项内力。用空间梁单元对结构进行离散分析时，计算结果直接给出截面的内力和变形。对于混凝土桥梁结构，理论计算和实验均证明，截面翘曲引起的正应力与按纯扭转理论所得应力值相比很小，通常不超过5%10%。一般按纯扭转理论进行分析可以满足设计要求。但对于钢箱梁，则必修考虑用翘曲扭转梁单元进行离散。当对宽箱梁桥分析时，空间梁单元法有很大的局限性，而且一维梁单元离散结构仅能得到系梁、拱、吊杆内力，而无法得到横梁的内力。板壳元。钢筋混凝土桥梁通常做成空间箱形结构，采用板、壳单元进行离散，当板壳单元相当细密时，可以包括桥梁结构的各种受力行为。板壳元法是分析桥梁上部结构最通用的方法，但在实际应用时，整理数据较为复杂工作量较大。因此在应用上受到很大限制。三维实体元。从整体结构中取出隔离体，按整体分析得到隔离体截面内力或位移条件作为隔离体边界条件，采用三维实体元对子结构进行分析。实体元可采用四面体、六面体以及等参元。三维实体元法在空间结构有限元分析中有一定的局限性，受计算机内存和计算速度的影响，实体单元离散整个结构构件难以实现，但是如果从整体结构中取出隔离体进行局部分析时，个隔离体的大小、内力、位移边界条件等因素对分析结果影响较大。梁格单元。梁格法时采用有限元技术分析桥梁上部结构的一种有效而又实用的方法。梁格法的特点是用一个等效梁格代替桥梁上部结构，分析此梁格的受力状态就可得实体桥梁的受力状态。该方法不仅适用于板式、梁板式及箱梁截面的上部结构，而且对分析弯、斜梁特别有效。二、如图示平面应力梁模型，每个节点2个自由度(u,v),试写出节点联系数组IT、单元联系数组LMT和对角元地址数组MAXA。2-1平面应力梁模型该平面应力梁模型单元信息数组：标识数组：节点联系数组: 三、求解线性方程组：,式中的为对称正定矩阵，其值为：矩阵为要求写出矩阵分解和回带求解过程。为了避免开方运算，利用改进的平方根法求解线性方程组。(1)矩阵分解其中L为单位下三角矩阵，D为对角矩阵。此时有：求解结果如下：(2)方程组的求解:则令则，求解结果如下：四、结构分析程序设计实践题：上机调试教材中的TRUSS程序，并要求给出具体算例。根据课本4246页，空间桁架的通用程序如下： PROGRAM TRUSS IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) CHARACTER NAME*40 COMMON /AT/A(18000000) COMMON /IAT/IA(2000000) A=0;IA=0 WRITE(*,*)INPUT FILE NAME? READ(*,*) NAME CALL OPENF(NAME) CALL DATAIN(NP,NE,NF,ND,NDF,NPF,NM,NR,NCF,& IME,INAE,IIT,ILMT,IMAXA,& IX,IY,IZ,IRR,IAE,IPF,ICKK) CALL FLMT(NP,NE,NN,NN1,NR,A(IRR+1),ND,NF,NDF,& IA(IME+1),IA(IIT+1),IA(ILMT+1) CALL FMAXA(NN1,NE,IA(ILMT+1),IA(IMAXA+1),NWK,NPF,NDF,& ICKK,IDIST,IFTOOL,IFF,IPP,ISG,ISM) CALL CONKB(NP,NE,NM,NWK,IA(IME+1),A(IX+1),& A(IY+1),A(IZ+1),A(IAE+1),IA(INAE+1),& IA(ILMT+1),IA(IMAXA+1),A(ICKK+1),NN1) CALL MKFORCE(NP,NF,NPF,NCF,NN,& IA(IIT+1),A(IPF+1),A(IPP+1),A(IFTOOL+1) CALL LDLT(A(ICKK+1),IA(IMAXA+1),NN,1,3,NWK,NN1) CALL RESOLVE(A(ICKK+1),A(IFTOOL+1),IA(IMAXA+1),NN,NWK,NN1) CALL DISPLS(NP,NE,NF,NPF,NM,NN,IA(IIT+1),A(IFTOOL+1),& A(IDIST+1),A(IAE+1),IA(IME+1),IA(INAE+1),A(IX+1),A(IY+1),& A(IZ+1), A(IPP+1),A(IFF+1),A(ISG+1),A(ISM+1) CALL DATAOUT(NP,NE,NPF,A(IDIST+1),A(IFF+1),A(ISG+1),A(ISM+1) CALL CLOSEF ENDC SUBROUTINE OPENF(NAME) CHARACTER NAME*40 NUM=0 DO I=1,40 IF(NAME(I:I).NE. )NUM=NUM+1 ENDDO OPEN(1,FILE=NAME(1:NUM), STATUS=UNKNOWN) OPEN(2,FILE=NAME(1:NUM)/.RES,STATUS=UNKNOWN) OPEN(3,FILE=NAME(1:NUM)/.ERO,STATUS=UNKNOWN) RETURN ENDCC SUBROUTINE CLOSEF CLOSE(1) CLOSE(2) CLOSE(3) RETURN ENDC SUBROUTINE DATAIN(NP,NE,NF,ND,NDF,NPF,NM,NR,NCF,& IME,INAE,IIT,ILMT,IMAXA,& IX,IY,IZ,IRR,IAE,IPF,ICKK) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) COMMON /AT/A(18000000) COMMON /IAT/IA(2000000) READ(1,*) NP,NE,NM,NR,NCF WRITE(2,701)NP,NE,NM,NR,NCF 701 FORMAT(/1X,#OUTPUT OF ORIGINAL INPUT INFORMATION#& /5X,Number of joints JOINTS=,I5& /5X,Number of elements ELEMENTS=,I5& /5X,Number of material property groups PROPERTY TYPES=,I5& /5X,Number of restrained joints RESTRAINTS=,I5& /5X,Number of concentrative forced joints NCF=,I5)C-FORM POINTER- NF=3 ND=2 NDF=ND*NF NPF=NP*NF IME=0 INAE=IME+2*NE IIT=INAE+NE ILMT=IIT+NF*NP IMAXA=ILMT+NDF*NE IX=0 IY=IX+NP IZ=IY+NP IRR=IZ+NP IAE=IRR+2*NR IPF=IAE+2*NM ICKK=IPF+4*NCF READ(1,*)(A(IX+I),A(IY+I),A(IZ+I),I=1,NP) WRITE(2,714)(I,A(IX+I),A(IY+I),A(IZ+I),I=1,NP) 714 FORMAT(/5X,GENERATED JOINT COORDINATES DATA& /1X, JOINT ,15X,X,13X,Y,13X,Z& /(4X,I5,3X,3(2X,E12.6) READ(1,*)(A(IAE+2*(I-1)+1),A(IAE+2*(I-1)+2),I=1,NM) READ(1,*)(IA(IME+2*(I-1)+1),IA(IME+2*(I-1)+2),IA(INAE+I),I=1,NE) WRITE(2,606)(I,A(IAE+2*(I-1)+1),A(IAE+2*(I-1)+2),I=1,NM) WRITE(2,607)(I,IA(IME+2*(I-1)+1),IA(IME+2*(I-1)+2),& IA(INAE+I),I=1,NE)606 FORMAT(/5X,ELEMENT MATERAIL PROPERTIES DATA& /2X,NO.,10X, E,10X,Ax& /(2X,I3,2(1X,E11.5)607 FORMAT(/5X,TRUSS ELEMENT DEFINITION DATA& /2X,NO.,10X,JOINT_1,10X,JOINT_2,10X,NAE& /(2X,I3,3(10X,I5) READ(1,*)(A(IRR+2*(I-1)+1), A(IRR+2*(I-1)+2),I=1,NR) WRITE(2,608)(A(IRR+2*(I-1)+1), A(IRR+2*(I-1)+2),I=1,NR)608 FORMAT(/5X,JOINT RESTRAINTS DATA& /2X, JOINT,10X,RESTRAINT,& /(2X,F7.0,10X,F9.3) READ(1,*)(A(IPF+4*(I-1)+J),J=1,4),I=1,NCF) WRITE(2,609)(A(IPF+4*(I-1)+J),J=1,4),I=1,NCF)609 FORMAT(/5X,CONCENTRATIVE FORCED JOINTS DATA& /2X, JOINT,10X,Fx, 10X,Fy,10X,Fz& /(2X,F7.0,3(1X,E12.6) RETURN ENDC SUBROUTINE MKFORCE(NP,NF,NPF,NCF,NN,IT,PF,PP,FTOOL) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION IT(NF,NP),PF(4,NCF),PP(NPF), FTOOL(NPF) PP=0;FTOOL=0 DO I=1,NCF NOD=PF(1,I) DO J=1,NF PP(NF*(NOD-1)+J)=PF(J+1,I) ENDDO ENDDO DO I=1,NP DO J=1,NF LAB=IT(J,I) IF(LAB.GT.0.AND.LAB.LE.NN) THEN FTOOL(LAB)=PP(NF*(I-1)+J) ENDIF ENDDO ENDDO RETURN ENDCC SUBROUTINE DATAOUT(NP,NE,NPF,DIST,FF,SG,SM) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION DIST(NPF),FF(NPF),SG(NE),SM(NE) WRITE(2,715)(I,(DIST(3*(I-1)+J),J=1,3),I=1,NP) 715 FORMAT(/5X,SOLVED JOINT DISPLACEMENTS DATA& /1X, JOINT ,3X,8X,Dx,12X,Dy,12X,Dz& /(4X,I5,3X,3(2X,E12.6) WRITE(2,716)(IE,SG(IE),SM(IE),IE=1,NE) 716 FORMAT(/5X,SOLVED ELEMENT INTERNAL FORCE DATA& /1X, ELEMENT ,3X,8X,Nx,8X,STRESS& /(4X,I5,3X,2(2X,F12.6) WRITE(2,717)(I,(FF(3*(I-1)+J),J=1,3),I=1,NP) 717 FORMAT(/5X,SOLVED JOINT REACTION DATA& /1X, JOINT ,3X,8X,Rx,12X,Ry,12X,Rz& /(4X,I5,3X,3(2X,f12.4) RETURN ENDCC SUBROUTINE FLMT(NP,NE,NN,NN1,NR,RR,ND,NF,NDF,ME,IT,LMT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N)C This program forms the joint&element numbering matrix IT&LMT DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR) NN=0;NN1=0;IT=0;LMT=0 N=0 DO I=1,NP C=0 DO K=1,NR KR=RR(1,K) IF(KR.EQ.I) C=RR(2,K) ENDDO NC=C C=C-NC DO J=1,NF C=C*10.0 L=C+0.1 C=C-L IF(L.EQ.0)THEN N=N+1 IT(J,I)=N ELSE IT(J,I)=0 ENDIF ENDDO ENDDO NN=N NN1=NN+1 DO IE=1,NE DO I=1,ND NI=ME(I,IE) DO J=1,NF LMT(I-1)*NF+J,IE)=IT(J,NI) ENDDO ENDDO ENDDO RETURN ENDC SUBROUTINE FMAXA(NN1,NE,LMT,MAXA,NWK,NPF,NDF,& ICKK,IDIST,IFTOOL,IFF,IPP,ISG,ISM)C This program forms the MDE address matrix MAXA of K IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION MAXA(NPF),LMT(NDF,NE) MAXA=0;NWK=0 MAXA(1)=1 DO I=2,NN1 IP=I-1 IG=IP DO IE=1,NE DO J=1,NDF IF(LMT(J,IE).EQ.IP) THEN DO K=1,NDF IF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE) ENDDO END IF ENDDO ENDDO MAXA(I)=MAXA(I-1)+IP-IG+1 ENDDO NWK=MAXA(NN1)-1 IDIST=ICKK+NWK IFTOOL=IDIST+NPF IFF=IFTOOL+NPF IPP=IFF+NPF ISG=IPP+NPF ISM=ISG+NE RETURN ENDC SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,& LMT,MAXA,CKK,NN1) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),& NAE(NE),LMT(6,NE),ME(2,NE),MAXA(NN1),& AKE(2,2),T(2,6),TT(6,2),AK(6,2),TAK(6,6) CKK=0 DO 10 IE=1,NE TAK=0 CALL FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE) CALL FT(IE,NP,NE,X,Y,Z,ME,T) CALL MAT(2,6,T,TT) AK=MATMUL(TT,AKE) TAK=MATMUL(AK,T) DO 220 I=1,6 DO 220 J=1,6 NI=LMT(I,IE) NJ=LMT(J,IE) IF(NJ-NI).GE.0.AND.NI*NJ.GT.0) THEN IJ=MAXA(NJ)+NJ-NI CKK(IJ)=CKK(IJ)+TAK(I,J) ENDIF 220 CONTINUE 10 CONTINUE RETURN ENDC SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NWK),MAXA(NNM) IF(NN.EQ.1) RETURN DO 200 N=1,NN KN=MAXA(N) KL=KN+1 KU=MAXA(N+1)-1 KH=KU-KL IF(KH)304,240,210210 K=N-KH IC=0 KLT=KU DO 260 J=1,KH KLT=KLT-1 IC=IC+1 KI=MAXA(K) ND=MAXA(K+1)-KI-1 IF(ND) 260,260,270270 KK=MIN0(IC,ND) C=0.0 DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L) A(KLT)=A(KLT)-C260 K=K+1240 K=N B=0.0 DO 300 KK=KL,KU K=K-1 KI=MAXA(K) C=A(KK)/A(KI) IF(ABS(C).LT.1.0E+07) GOTO 290 WRITE(IOUT,2010) N,C STOP290 B=B+C*A(KK)300 A(KK)=C A(KN)=A(KN)-B304 IF(A(KN) 310,310,200310 IF(ISH.EQ.0) GOTO 320 IF(A(KN).EQ.0.0) A(KN)=-1.0E-16 GOTO 200320 WRITE(IOUT,2000) N,A(KN) STOP200 CONTINUE RETURN2000 FORMAT(/ Stop-stiffness matrix not positive + definite,/,nonpositive pivot for equation, + I4,/, pivot =,E20.10)2010 FORMAT(/, Stop-sturm sequence check failed + because of multiplier growth for column + number,I4,/, Multiplier = ,E20.8) ENDCC SUBROUTINE RESOLVE(A,V,MAXA,NN,NWK,NNM) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NWK),V(NN,1),MAXA(NNM) NIP=1 DO IP=1,NIP DO 400 N=1,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF(KU-KL) 400,410,410410 K=N C=0.0 DO 420 KK=KL,KU K=K-1420 C=C+A(KK)*V(K,IP) V(N,IP)=V(N,IP)-C400 CONTINUE DO 480 N=1,NN K=MAXA(N)480 V(N,IP)=V(N,IP)/A(K) IF(NN.EQ.1)RETURN N=NN DO 500 L=2,NN KL=MAXA(N)+1 KU=MAXA(N+1)-1 IF(KU-KL) 500,510,510510 K=N DO 520 KK=KL,KU K=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1 ENDDO RETURN ENDCCcc SUBROUTINE DISPLS(NP,NE,NF,NPF,NM,NN,IT,FTOOL,& DIST,AE,ME,NAE,X,Y,Z,PP,FF,SG,SM) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION IT(NF,NP),DIST(NPF),FTOOL(NPF),T(2,6),& TT(6,2),AE(2,NM),ME(2,NE),NAE(NE),UE(6),U(2), & AKE(2,2),FE1(2),FE(6),FF(NPF),X(NP),Y(NP),Z(NP),& PP(NPF),SG(NE),SM(NE) SG=0;SM=0;FF=0 DO I=1,NP DO J=1,NF LAB=IT(J,I) IF(LAB.EQ.0) THEN DIST(NF*(I-1)+J)=0.0 ELSEIF(LAB.GT.0.AND.LAB.LE.NN) THEN DIST(NF*(I-1)+J)=FTOOL(LAB) ENDIF ENDDO ENDDO DO IE=1,NE N1=ME(1,IE);N2=ME(2,IE) UE=0 DO J=1,NF UE(J)=DIST(NF*(N1-1)+J) UE(NF+J)=DIST(NF*(N2-1)+J) ENDDO CALL FT(IE,NP,NE,X,Y,Z,ME,T) CALL FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE) U=MATMUL(T,UE) FE1=MATMUL(AKE,U) CALL MAT(2,6,T,TT) FE=MATMUL(TT,FE1) DO J=1,NF FF(NF*(N1-1)+J)=FF(NF*(N1-1)+J)+FE(J) FF(NF*(N2-1)+J)=FF(NF*(N2-1)+J)+FE(NF+J) ENDDO ISW=NAE(IE) AO=AE(2,ISW) SG(IE)=FE1(2) SM(IE)=FE1(2)/AO ENDDO DO I=1,NPF FF(I)=FF(I)-PP(I)ENDDO RETURN ENDC SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),NAE(NE),AE(2,NM),AKE(2,2) N1=ME(1,IE) N2=ME(2,IE) X1=X(N1);Y1=Y(N1);Z1=Z(N1) X2=X(N2);Y2=Y(N2);Z2=Z(N2) BL=SQRT(X2-X1)*2+(Y2-Y1)*2+(Z2-Z1)*2) NMI=NAE(IE) E0=AE(1,NMI);A0=AE(2,NMI) C=E0*A0/BL AKE(1,1)=C AKE(1,2)=-C AKE(2,1)=-C AKE(2,2)=C RETURN ENDC SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6) T=0 N1=ME(1,IE);N2=ME(2,IE) X1=X(N1);Y1=Y(N1);Z1=Z(N1) X2=X(N2);Y2=Y(N2);Z2=Z(N2) BL=SQRT(X2-X1)*2+(Y2-Y1)*2+(Z2-Z1)*2) CX=(X2-X1)/BL CY=(Y2-Y1)/BL CZ=(Z2-Z1)/BL T(1,1)=CX;T(2,4)=CX T(1,2)=CY;T(2,5)=CY T(1,3)=CZ;T(2,6)=CZ RETURN END C SUBROUTINE MAT(M,N,A,B) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION A(M,N),B(N,M) DO I=1,M DO J=1,N B(J,I)=A(I,J) END DO END DO RETURN END C该程序在Microsoft Developer Studio上的Fortran PowerStation 4.0 上调试，该程序运行无误。表4-1 程序说明输入数据控制数据NP结构离散节点总数NE结构离散单元总数NM结构单元不同特征数类总数NR结构受约束节点总数NCF结构外荷载作用的总数节点坐标A(IX+I)节点X坐标A(IY+I)节点Y坐标A(IX+I)节点Z坐标材料特性A(IAE+2*(I-1)+1)，A(IAE+2*(I-1)+2)单元信息IA(IME+2*(I-1)+1)，IA(IME+2*(I-1)+2)单元信息ME约束信息A(IRR+2*(I-1)+1)， A(IRR+2*(I-1)+2)约束信息RR荷载信息A(IPF+4*(I-1)+J)外荷载数据PF输出数据节点位移DIST(3*(I-1)+J)单元内力SG(IE)单元内力SM(IE)单元截面应力节点约束力FF(6*(I-1)+J)算例一：现将教材Page20 图2-7 空间桁架结构作为算例进行计算。如图所示：图4-1空间桁架结构图说明：(1) 所有材料的弹性模量为2.20108kPa，截面积为0.008m2；(2) 桁架尺寸图中已标出单位:m；(3) 在节点7和节点8分别施加三个力，方向分别沿X轴正向、Y轴正向、Z轴负向，大小均为10kN。输入文件信息如下：控制数据：815142节点坐标：000400800040440880402442材料特征：2.20E+080.008单元信息：121231141251361451561171271371481581681781281约束信息：10.11130.11140.11160.111荷载信息：71010-1081010-10说明：叙述说明语句不得输入文件。计算结果在工程文件夹的.RES文件中，结果如下： #OUTPUT OF ORIGINAL INPUT INFORMATION# Number of joints JOINTS= 8 Number of elements ELEMENTS= 15 Number of material property groups PROPERTY TYPES= 1 Number of restrained joints RESTRAINTS= 4 Number of concentrative forced joints NCF= 2 GENERATED JOINT COORDINATES DATA JOINT X Y Z 1 .000000E+00 .000000E+00 .000000E+00 2 .400000E+01 .000000E+00 .000000E+00 3 .800000E+01 .000000E+00 .000000E+00 4 .000000E+00 .400000E+01 .000000E+00 5 .400000E+01 .400000E+01 .000000E+00 6 .800000E+01 .800000E+01 .000000E+00 7 .400000E+01 .000000E+00 .200000E+01 8 .400000E+01 .400000E+01 .200000E+01 ELEMENT MATERAIL PROPERTIES DATA NO. E Ax 1 .22000E+09 .80000E-02 TRUSS ELEMENT DEFINITION DATA NO. JOINT_1 JOINT_2 NAE 1 1 2 1 2 2 3 1 3 1 4 1 4 2 5 1 5 3 6 1 6 4 5 1 7 5 6 1 8 1 7 1 9 2 7 1 10 3 7 1 11 4 8 1 12 5 8 1 13 6 8 1 14 7 8 1 15 2 8 1 JO