结构程序设计

1、 和？7夫孝 CENTRAL SOUTH UNIVERSITY 结构分析程序设计基础 姓 名： 学 号： 专 业： 授课教师： 指导教师： 2016年5月20日 一、结合自身所学专业，是阐述利用有限单元法进行结构分析时所常用离散模 型及其特点。 应用有限单元法求解土木工程问题时常用的离散模型有空间梁单元、板壳 元、三维实体元及梁格单元。 空间梁单元。空间梁单元是一种常用的单元，除了桥梁上部结构外，框架系统 和动力联动装置都可以采用梁单元模拟。对于一个采用梁单元模拟的构件，一般 它在一个方向上的尺寸最少是另两个方向上的10倍。梁单元除了承担拉力和压 力，还可以承担剪力和弯矩。三维梁单元每个节点有

2、 6个自由度，即沿轴的位移 和绕这3个轴的转角。根据结构受载后截面是否保持平截面， 可分为自由扭转和 翘曲扭转。后者考虑了受载后横截面双力矩阵和翘曲扭矩两项内力。用空间梁单 元对结构进行离散分析时，计算结果直接给出截面的内力和变形。 对于混凝土桥 梁结构，理论计算和实验均证明，截面翘曲引起的正应力与按纯扭转理论所得应 力值相比很小，通常不超过5%10%。一般按纯扭转理论进行分析可以满足设计 要求。但对于钢箱梁，则必修考虑用翘曲扭转梁单元进行离散。当对宽箱梁桥分 析时，空间梁单元法有很大的局限性，而且一维梁单元离散结构仅能得到系梁、 拱、吊杆内力，而无法得到横梁的内力。 板壳元。钢筋混凝土桥梁通

3、常做成空间箱形结构， 采用板、壳单元进行离散， 当板壳单元相当细密时，可以包括桥梁结构的各种受力行为。 板壳元法是分析桥 梁上部结构最通用的方法，但在实际应用时，整理数据较为复杂工作量较大。因 此在应用上受到很大限制。 三维实体元。从整体结构中取出隔离体，按整体分析得到隔离体截面内力或 位移条件作为隔离体边界条件，采用三维实体元对子结构进行分析。实体元可采 用四面体、六面体以及等参元。三维实体元法在空间结构有限元分析中有一定的 局限性，受计算机内存和计算速度的影响，实体单元离散整个结构构件难以实现， 但是如果从整体结构中取出隔离体进行局部分析时，个隔离体的大小、内力、位 移边界条件等因素对分析

4、结果影响较大。 梁格单元。梁格法时采用有限元技术分析桥梁上部结构的一种有效而又实用 的方法。梁格法的特点是用一个等效梁格代替桥梁上部结构，分析此梁格的受力 状态就可得实体桥梁的受力状态。 该方法不仅适用于板式、梁板式及箱梁截面的 上部结构，而且对分析弯、斜梁特别有效。 二、如图示平面应力梁模型，每个节点 2个自由度（u,v）,试写出节点联系数组 IT、单元联系数组LMT和对角元地址数组MAXA。 2-1平面应力梁模型 ME（3,24） 112 2 6 2 7 3 标识数组： 4 10 11 12 8 14 10 15 11 16 12 18 14 19 15 20 16 该平面应力梁模型单元信

5、息数组: 3 55667799 10 10 11 11 13 13 14 14 15 15 8 910 10 11 11 12 13 14 14 15 15 16 17 18 18 19 19 20 1 ID（2,0）1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 节点联系数组： 013579 IT（2,0）0246810 11 13 15 17 12 14 16 18 19 21 23 25 27 29 20 22 24 26 28 30 0 0 7 8 9 10 1 1 1 2 0 0 9 10 1 2 1 4 1 2 9 10 11 12 1 7 1 2 11 12

6、3 4 3 11 3 4 11 12 13 14 3 14 7 18 3 4 13 14 5 6 7 19 7 8 15 16 17 18 7 21 7 8 17 18 9 10 7 24 9 10 17 18 19 20 1 28 1 39 9 10 19 20 11 12 3 51 11 12 19 20 21 22 3 62 11 12 21 22 13 14 7 74 LMT (6,24)丁 mi 7 MAXA 83 15 16 23 24 25 26 7 93 15 16 25 26 17 18 7 104 17 18 25 26 27 28 9 116 17 18 27 28 19

7、 20 9 127 11 139 19 20 27 28 29 30 11 150 1920 2930212215162 23 24 0 0 0 0 15 171 23 24 0 0 25 26 15 181 15 192 25 26 0 0 0 0 17 204 25 26 0 0 27 28 17 215 27 28 0 0 0 0 19 227 19 238 27 28 0 0 29 30 19 250 、求解线性方程组：A X B ,式中的A为对称正定矩阵，其值为: 45 2 -13 0 0 0 0 53 0 0 0 0 0 102 51 0 -17 0 A 84 0 0 0 对称 6

8、 0 0 31 0 29 B矩阵为 要求写出矩阵分解和回带求解过程。 100 200300 400 500 500 2000 为了避免开方运算，利用改进的平方根法求解线性方程组。 (1)矩阵分解 A LDLt 其中 1l为单位下三角矩阵, d为对角矩阵。 此时有： dj j 1 ajj1 jkd k 1 k lij j 1 . (aijhkl jkdk) d j, i k 1/ j 1, j 2, n. 求解结果如下： 1.0000 0 0 0 0 0 0 0.0444 1.0000 0 0 0 0 0 0.2889 0.0109 1.0000 0 0 0 0 L 0 0 0.5191 1.0

9、000 0 0 0 0 0 0 0 1.0000 0 0 0 0 0.1730 0.1534 0 1.0000 0 0 0 0 0 0 0 1.0000 45.0000 0 0 0 0 0 0 0 52.9111 0 0 0 0 0 0 0 98.2381 0 0 0 0 D 0 0 0 57.5235 0 0 0 0 0 0 0 6.0000 0 0 0 0 0 0 0 26.7041 0 0 0 0 0 0 0 29.0000 (2)方程组的求解： A XB 则LDLt XB 令DLt X Y 则L YB，Lt X Dt 丫 k 1 ykbklkjyj,k 1,2,n. j 1 n Xky

10、k.dk lkjyj,k n,n 1,1. j k 1 求解结果如下: 100.0 195.6 326.8 230.4 500.0 521.2 2000.0 X 3.84603.62856.1792 1.0103 83.3333 19.5176 68.9655 四、结构分析程序设计实践题：上机调试教材中的TRUSS程序，并要求给出 具体算例。 根据课本4246页，空间桁架的通用程序如下： PROGRAM TRUSS IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) CHARACTER NAME*40 COMMON /A T/A(1800000

11、0) 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, 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.O.AND.LAB .L E.NN) THEN FTOOL(LAB)=PP(NF

12、*(I-1)+J) ENDIF ENDDO ENDDO RETURN END C C 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)(l,(DIST(3*(l-1)+J),J=1,3),l=1,NP) 715 FORMA T(/5X,SOLVED JOINT DISPLACEMENTS DATA NN1=0;IT=0;LMT=0 N=0 DO I

13、=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.O)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 END C SUBROUTINE FMAXA(

14、NN1,NE,LMT,MAXA,NWK,NPF,NDF, NWK=0 MAXA(1)=1 DO l=2,NN1 IP=I-1 IG=IP DO IE=1,NE DO J=1,NDF IF(LMT(J,IE).EQP) THEN DO K=1,NDF IF(LMT(K,IE).GT.O.AND 丄 MT(K,IE ).L E.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

15、IPP=IFF+NPF ISG=IPP+NPF ISM=ISG+NE RETURN END C SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y ,Z,AE,NAE, 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 2= ME(1,IE);N2=ME(2,IE) UE=0 D

16、O 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=MA TMUL(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) E

17、NDDO 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 END C 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,I

18、E) X仁 X(N1);Y 仁丫(N1);Z 仁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 END C SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION X(N

19、P),Y(NP),Z(NP),ME(2,NE),T(2,6) T=0 2= ME(1,IE);N2=ME(2,IE) X1=X(N1);Y 1=Y(N1);Z 1=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

20、 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) 节点

21、X坐标 A(IY+I) 节点Y坐标 A(IX+I) 节点Z坐标 材料特性 A(IAE+2*(l-1)+1)，A(IAE+2*(l-1)+2) 单元信息 IA(IME+2*(l-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*(l-1)+J) 说明： 所有材料的弹性模量为2.20X 108kPa,截

22、面积为0.008m2; (2)桁架尺寸图中已标出单位：m； 在节点7和节点8分别施加三个力，方向分别沿 X轴正向、丫轴正向、Z轴 负向，大小均为10kN。 输入文件信息如下： 控制数据： 8151 42 节点坐标： 000 4 00 800 040 4 40 880 402 442 材料特征： 2.20E+080.008 单元信息： 1 2 1 2 3 1 1 41 2 51 3 61 4 51 5 61 1 71 2 71 3 71 4 81 5 81 6 8 1 7 81 28 1 约束信息： 10.111 3 0.111 4 0.111 6 0.111 荷载信息： 7 1010-10 8

23、 10 10 -10 说明：叙述说明语句不得输入文件。 计算结果在工程文件夹的.RES文件中，结果如下: #OUTPUT OF ORIGINAL INPUT INFORMATION# JOINTS= ELEMENTS= PROPERTY TYPES= 1 RESTRAINTS= 4 NCF= 2 8 15 Number of joints Number of eleme nts Number of material property groups Number of restrai ned joi nts Number of concen trative forced joints GENER

24、ATED 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 .20000

25、0E+01 8 .400000E+01 .400000E+01 .200000E+01 ELEMENT MA TERAIL PROPERTIES DA TA NO.EAx 1 .22000E+09.80000E-02 TRUSS ELEMENT DEFINITION DATA NO.JOINT_1JOINT_2NAE 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 JOINT RESTRAIN

26、TS DA TA JOINT RESTRAINT 1. .111 3. .111 4. .111 6. .111 CONCENTRA TIVE FORCED JOINTS DATA JOINT FxFyFz 7. .100000E+02.100000E+02 -.100000E+02 8. .100000E+02.100000E+02 -.100000E+02 SOLVED JOINT DISPLACEMENTS DATA JOINT Dx Dy Dz 1 .000000E+00 .000000E+00 .000000E+00 2 .000000E+00 .365790E-04 -.51043

27、3E-04 3 .000000E+00 .000000E+00 .000000E+00 4 .000000E+00 .000000E+00 .000000E+00 5 -.757576E-05 .290032E-04 -.849140E-04 6 .000000E+00 .000000E+00 .000000E+00 158812E-04.868291E-04-.529372E-04 .212821E-04.641018E-04-.849140E-04 SOLVED ELEMENT INTERNAL FORCE DA TA ELEMENT Nx STRESS 1 .000000 .000000

28、 2 .000000 .000000 3 .000000 .000000 4 -3.333333 -416.666667 5 .000000 .000000 6 -3.333333 -416.666667 7 -4.714045 -589.255651 8 -3.726780 -465.847495 9 -1.666667 -208.333333 10 -14.907120 -1863.389981 11 -7.453560 -931.694991 12 .000000 .000000 13 -25.000000 -3125.000000 14 -10.000000 -1250.000000

29、15 3.726780 465.847495 SOLVED JOINT REACTION DA TA JOINT Rx Ry Rz 1 3.3333 .0000 1.6667 2 .0000 .0000 .0000 3 -13.3333 .0000 6.6667 4 10.0000 .0000 3.3333 5 .0000 .0000 .0000 6 -20.0000 -20.0000 8.3333 7 .0000 .0000 .0000 8 .0000 .0000 .0000 整理如下： 表4-2节点位移表 JOINT Dx (mm) Dy (mm) Dz (mm) 1 0 0 0 2 0

30、0.036579 -0.051043 3 0 0 0 4 0 0 0 5 -0.00757576 0.0290032 -0.084914 6 0 0 0 7 0.0158812 0.0868291 -0.052937 8 0.0212821 0.0641018 -0.084914 ELEMENT Nx (kN) STRESS (kPa) 1 0 0 2 0 0 3 0 0 4 -3.333333 -416.6667 5 0 0 6 -3.333333 -416.6667 7 -4.714045 -589.2557 8 -3.72678 -465.8475 9 -1.666667 -208.33

31、33 10 -14.90712 -1863.39 11 -7.45356 -931.695 12 0 0 13 -25 -3125 14 -10 -1250 15 3.72678 465.8475 表4-3单元计算结果 JOINT Rx (kN) Ry (kN) Rz (kN) 1 3.3333 0 1.6667 2 0 0 0 3 -13.3333 0 6.6667 4 10 0 3.3333 5 0 0 0 6 -20 -20 8.3333 7 0 0 0 8 0 0 0 表4-4支座反力 算例二： 图示空间桁架为圆形桁架穹顶，几何尺寸，坐标系位置，及荷载作用位置如 图所示。 说明： 所有

32、材料的弹性模量为2.20X 108kPa,截面积为0.06m2; (2)桁架尺寸单位：m； 在节点9、节点12和节点13分别施加一个力，方向分别沿 X轴负向、X轴 正向、丫轴负向，分别为10kN、10kN、50kN。 图4-2桁架立面图 4 图4-3桁架平面图 输入文件信息如下： 控制数据： 1324163 节点坐标: 0-10 50 43.3 -10 25 43.3 -10 -25 0 -10 -50 -43.3 -10 -25 -43.3 -10 25 -12.5 -2 21.65 12.5 -2 21.65 25-2 0 12.5 -2 -21.65 -12.5 -2 -21.65 -2

33、5 -2 0 0 0 0 材料特性： 2.2E+08 0.06 单元信息: 2 9 1 3 9 1 3 10 1 4 10 1 4 11 1 5 11 1 5 12 1 6 12 1 6 7 1 7 8 1 8 9 1 9 10 1 10 11 1 11 12 1 7 12 1 7 13 1 8 13 1 9 13 1 10 13 1 11 13 1 12 13 1 约束信息： 10.111 2 0.111 3 0.111 4 0.111 5 0.111 6 0.111 荷载信息： 9-10 00 1210 00 130-500 计算结果整理如下: 4-5节点位移表 JOINT Dx(m) Dy(m) Dz(m) 1 0.00E+00 0.00E+00 0.00E+00 2 0.00E+00 0.00E+00 0.00E+00 3 0.00E+00 0.00E+00 0.00E+00 4 0.00E+00 0.00E+00 0.00E+00 5 0.00E+00 0.00E+00 0.00E+00 6 0.00E+00 0