空间桁架结构程序设计(Fortran)学习资料
第二章桁架结构的有限元
2-8 计算杆件内力
计算出单元节点位移{ 计算出单元节点位移 ui,vi,uj,vj}T，可计算出单元两端的 节点力和内力。 节点力和内力。 轴向力： 轴向力：
1 Ui 0 0 EA = lu i 0 0 0 v i u j 0 1 0 0 0 0 v j
1、建立如图所示的杆系结构； 、建立如图所示的杆系结构； 2、定义单元类型：LINK1 、定义单元类型： 3、定义材料弹性模量EX 、定义材料弹性模量 4、定义实常数：杆的截面积0.01 、定义实常数：杆的截面积 5、划分网格：一个杆为一个单元 、划分网格： 6、定义约束 、 7、施加载荷 、 8、进行求解 、 9、观察变形图、列出节点位移值 、观察变形图、
δ （-4,3） 2 ① 3 （0,5） 4 （4,3）
②
③ EA=5e6N 1 P=1000N x
采用ANSYS 分析，计算节点的位移、反作用力和桁架系 统的应力。 几何参数及载荷如图3-10所示，杆的弹性模量E 为200Gpa， 横截面面积A 为3250mm2。
图3-10桥梁桁架模型
2-7添加约束 求解方程 添加约束
约束条件2：节点 水平位移为 水平位移为u 约束条件 ：节点n水平位移为 n=un*≠0： ： 在整体刚度矩阵K中 在整体刚度矩阵 中，与un相对应的行与列中主对角线元 乘以一个大数A，在右边向量F中 素K2n-1,2n-1乘以一个大数 ，在右边向量 中，与un相对应 的行元素改为AK2n-1,2n-1 un*，其他元素不变 ； 的行元素改为 经过这样修改后的位移法基本方程 K *δ = F * 可解出节 点位移δ 点位移
2-2
坐标转换的概念
在用有限元法计算中，第一步是将结构离散 在用有限元法计算中， 将结构离散成有限个单元， 化，将结构离散成有限个单元，一般一个杆 作为一个单元；在该单元的坐标系（ 作为一个单元；在该单元的坐标系（局部坐 标系）中建立单元刚度矩阵， 标系）中建立单元刚度矩阵，所有的单元刚 度矩阵（局部坐标系下） 度矩阵（局部坐标系下）需要整和成总体刚 度矩阵（整体坐标系下），即每个单元对整 度矩阵（整体坐标系下），即每个单元对整 ）， 体的贡献， 体的贡献，在整和过程中需要根据局部坐标 系与整体坐标系之间的关系（ 系与整体坐标系之间的关系（称为坐标转换 矩阵）进行坐标转换。 矩阵）进行坐标转换。
桁架结构设计
FN2 =3.33FP
FN3 =-0.50FP
截面单杆 截面法取出的隔离体，不管其上 有几个轴力，如果某杆的轴力可以通过列一 个平衡方程求得，则此杆称为截面单杆。 可能的截面单杆通常有相交型和平行型两种 形式。
小结： 熟练掌握 计算桁架内力的基 本方法： 结点法和截面法 采取最简捷的途径计算桁架内力 能够分析和计算组合结构的内力 尤其注意区分二力杆和非二力杆
2.5.2 桁架结构的分类：
一、根据维数分类 1. 平面（二维）桁架（plane truss） ——所有组成桁架的杆件以及荷载的作 用线都在同一平面内
2. 空间（三维）桁架（space truss） ——组成桁架的杆件不都在同一平面内
二、按外型分类 1. 平行弦桁架 2. 三角形桁架
3. 抛物线桁架
FN=0
FN1=0
判断结构中的零杆
FP FP FP/2
FP/ 2
FP
2.5.3
结点法（nodal analysis method）
以只有一个结点的隔离体为研究对象，用 汇交力系的平衡方程求解各杆的内力的方法 例1. 求以下桁架各杆的内力
0
-33 34.8 19 19
YNAD CD 0.5 X NAD AC 1.5
Y
S
α
X L Ly
α
S
Lx
S X Y L Lx Ly
2. 结点单杆 以结点为平衡对象能仅用一个方程 求出内力的杆件，称为结点单杆（nodal single bar）。 利用这个概念，根据荷载状况可判断此杆内力是 否为零。 3. 零杆 零内力杆简称零杆（zero bar）。
FN2=0
FN=0
0
-33 34.8 19 -8
桁架有限元程序流程（有限元课程设计）
有限单元法课程设计有限单元法是基于连续介质力学基础上发展起来的,目前使用最广泛的数值计算方法。
有限单元法解决问题的前提是各单元相邻边界的位移协调。
有限单元法解决问题的前提是各单元相邻边界的位移协调。
有限单元有限单元法将连续的求解域离散为一组有限个单元组成的组合体,由细分单元去逼近求解域,由于单元的不同连接方式和形式各异的单元形状由于单元的不同连接方式和形式各异的单元形状,,因此可以适应几何形状复杂的求解区域杂的求解区域;;第二第二,,利用每一个单元内的近似函数利用每一个单元内的近似函数((形函数形函数))来表示全求解域上待求的未知场函数待求的未知场函数,,把一个连续的无限自由度问题变成离散的有限自由度问题,只要求出单元结点的物理量只要求出单元结点的物理量,,就可以确定单元组合体上的其他未知场函数就可以确定单元组合体上的其他未知场函数,,如果选择合适的形函数选择合适的形函数,,随着网格密度的减小随着网格密度的减小,,近似解将逐步趋向精确解近似解将逐步趋向精确解;;第三第三,,有限单元法计算得到的总体刚度矩阵为稀疏带状矩阵,这样借助于电子计算机存储和计算的效率大大提高计算的效率大大提高,,便于处理大规模问题。
便于处理大规模问题。
从上述有限单元法的特性可知从上述有限单元法的特性可知,,其计算原理简单其计算原理简单,,但由于单元连接方式和单元形状的多元化元形状的多元化,,以及近似函数的选择合适与否以及近似函数的选择合适与否,,使得有限元法在针对具体问题求解时比较烦琐求解时比较烦琐,,正是基于这样的应用背景正是基于这样的应用背景,,本论文提出了一种更简单实用的单元模型—平面等效桁架单元模型。
最后最后,,编制有限元分析程序编制有限元分析程序,,将这种桁架单元模型运用于钢筋混凝土结构中型运用于钢筋混凝土结构中,,模型中混凝土采用等效桁架单元模型中混凝土采用等效桁架单元,,钢筋采用一维杆单元单元,,利用混凝土等效的应力应变关系对各种构件进行弹塑性分析,并试探性的提出了单元破坏准则。
有限元方法与MATLAB程序设计 第3章 桁架和刚架
Fe keUe
Ue TUe
结构坐标系下的 单元刚度方程
9
结构坐标系下单元刚度矩阵 ke T TkeT
c cos s sin
y
Fe keUe
y
i (xi , yi )
(xj, yj ) x j
x
结构坐标系
cos sin 0
0 1 0 1 0 e cos sin 0
0
ke
EA sin
Fyi
Fxj Fyj
EA l
0
1
0
0 0 0
0 1 0
0
vi
0 0
u v
j j
单元坐标系下的 1 0 1 0
单元刚度矩阵
ke
EA l
0
1
0 0
0 1
0 0
0
0
0
0
7
§3.1.3 整体坐标系下的单元刚度方程
Fxi Fxi cos Fyi sin Fyi Fxi sin Fyi cos Fxj Fxj cos Fyj sin
是结构坐标系x轴正方 向至单元坐标轴x 的角度
y
x
y Fxi i
Fxi Fyi Fyi
x
结构（整体）坐标系
Fyj Fxj sin Fyj cos
Fxi Fyi
Fxj
Fyj
e
cos
sin
0
0
sin cos
0 0
0 0
cos sin
0
e
Fxi
e0Fyi源自sin cosFxj
Fyj
Fe TFe Ue TUe
cos
T sin
0
0
sin cos
第二章桁架结构的有限元
2-2
坐标转换的概念
在用有限元法计算中，第一步是将结构离散 在用有限元法计算中， 将结构离散成有限个单元， 化，将结构离散成有限个单元，一般一个杆 作为一个单元；在该单元的坐标系（ 作为一个单元；在该单元的坐标系（局部坐 标系）中建立单元刚度矩阵， 标系）中建立单元刚度矩阵，所有的单元刚 度矩阵（局部坐标系下） 度矩阵（局部坐标系下）需要整和成总体刚 度矩阵（整体坐标系下），即每个单元对整 度矩阵（整体坐标系下），即每个单元对整 ）， 体的贡献， 体的贡献，在整和过程中需要根据局部坐标 系与整体坐标系之间的关系（ 系与整体坐标系之间的关系（称为坐标转换 矩阵）进行坐标转换。 矩阵）进行坐标转换。
y
y R
x
θ O
u
v
u = u cos θ + v sin θ v = −u sin θ + v cos θ
角度θ： 轴逆时针转到 轴为正值 角度 ： x轴逆时针转到 x
v
u
x
2-4 整体坐标系下的单元刚度矩阵 整体坐标系下的单元刚度矩阵
写成矩阵的形式为： 写成矩阵的形式为：
u cos θ = v − sin θ sin θ u v cos θ
有限元基础与ANSYS入门 有限元基础与ANSYS入门 ANSYS
Finite Element Foundation and ANSYS introduction
机械工程系
第二章 桁架结构有限元
第二章 桁架结构有限元的步骤
桁架结构是指结构由许多细长杆件构成的结构系统， 桁架结构是指结构由许多细长杆件构成的结构系统，且 杆件的弯曲刚度小， 杆件的弯曲刚度小，杆件的变形主要是轴向变形
(完整word版)空间桁架结构程序设计(Fortran)
空间桁架静力分析程序及算例1、变量及数组说明2、空间桁架结构有限元分析程序源代码!主程序（读入文件，调用总计算程序，输出结果）CHARACTER IDFUT*20,OUTFUT*20WRITE(*,*) 'Input Data File name:'READ (*,*)IDFUTOPEN (11,FILE=IDFUT,STATUS='OLD')WRITE(*,*) 'Output File name:'READ (*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')WRITE(12,*)'*****************************************'WRITE(12,*)'* Program for Analysis of Space Trusses *'WRITE(12,*)'* School of Civil Engineering CSU *'WRITE(12,*)'* 2012.6.25 Designed By MuZhaoxiang *'WRITE(12,*)'*****************************************'WRITE(12,*)' 'WRITE(12,*)'*****************************************'WRITE(12,*)'*************The Input Data****************'WRITE(12,*)'*****************************************'WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100 FORMAT(6X,'The General Information'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',& 5X,'NCF',5X,'ND')110 FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!********************************************************************!总计算程序SUBROUTINE ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR), NAE(NE),&AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),&PP(NPF),FF(NPF),SG(NE),SM(NE)READ(11,*)(X(I),Y(I),Z(I),I=1,NP)READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)READ(11,*)(RR(1,J),RR(2,J),J=1,NR)READ(11,*)(AE(1,J),J=1,2)WRITE(12,120)WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)WRITE(12,130)WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)WRITE(12,140)WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)WRITE(12,150)WRITE(12,151)(AE(1,J),J=1,2)IF(NCF/=0)THENREAD(11,*)((PF(I,J),I=1,4),J=1,NCF)WRITE(12,160)WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)ENDIF120 FORMAT(/6X,'The Information of Joints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')121 FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,'The Information of Members'/2x,'Member',2X,'START',4X,'END',6X,'NAE')131 FORMAT(1X,I4,3I8)140 FORMAT(/6X,'The Information of SUPPORTS'/2x,'Joint',5X,'S')141 FORMAT(1X,I4,F8.3)150 FORMAT(/6X,'The Information of Sections'/4x,'E0',8X,'A0')151 FORMAT(1X,1PE8.2,F8.4)160 FORMAT(/6X,'The Loading at Joints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')161 FORMAT(1X,I4,3F8.2)CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)ISH=1CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALL REBACK(CKK,V,MAXA,NN,NWK,NNM)CALL DISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF)END!********************************************************************!矩阵转置子程序SUBROUTINE MAT(M,N,A,B)DIMENSION A(M,N),B(N,M)DO I=1,MDO J=1,NB(J,I)=A(I,J)END DOEND DORETURNEND!单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)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/BLAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND!单元坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)T=0N1=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)/BLCY=(Y2-Y1)/BLCZ=(Z2-Z1)/BLT(1,1)=CX;T(2,4)=CXT(1,2)=CY;T(2,5)=CYT(1,3)=CZ;T(2,6)=CZRETURNEND!生成单元联系数组LMTSUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0DO I=1,NPC=0DO K=1,NRKR=RR(1,K)IF(KR.EQ.I) C=RR(2,K)ENDDONC=C !NC=0,提取了整数部分C=C-NC !C=0.***,例如C=0.111DO J=1,NFC=C*10.0 !例如C=1.21L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位 !上的数字,这里"+0.1"是为了防止四舍五入是出现错误 C=C-LIF(L.EQ.0)THENN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DO IE=1,NEDO I=1,NDNI=ME(I,IE)DO J=1,NFLMT((I-1)*NF+J,IE)=IT(J,NI)ENDDOENDDOENDDORETURNEND!二维总刚中对角线元地址数组SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)DIMENSION MAXA(NPF),LMT(NDF,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NNMIP=I-1IG=IPDO IE=1,NEDO J=1,NDFIF(LMT(J,IE).EQ.IP) THENDO K=1,NDFIF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE)ENDDOEND IFENDDOENDDOMAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK= MAXA(NNM)-1RETURNEND!生成一维存储结构总刚度矩阵SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),& MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)CKK=0DO 10 IE=1,NETAK=0CALL 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,6DO 220 J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IF((NJ-NI).GE.0.AND.NI*NJ.GT.0) THENIJ=MAXA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIF220 CONTINUE10 CONTINUERETURNEND!生成荷载矩阵SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF)V=0PP=0DO I=1,NFDO J=1,NPDO K=1,NCFIF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J))=PF(I+1,K)ENDIFENDDOENDDOENDDODO K=1,NCFDO I=1,NPIF(I.EQ.PF(1,K))THENPP(NF*(I-1)+1)=PF(2,K)PP(NF*(I-1)+2)=PF(3,K)PP(NF*(I-1)+3)=PF(4,K)ENDIFENDDOENDDORETURNEND!对一维结构总刚度矩阵进行矩阵分解（LDLT）SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM) DIMENSION A(NWK),MAXA(NNM)IF(NN.EQ.1) RETURNDO 200 N=1,NNKN=MAXA(N)KL=KN+1KU=MAXA(N+1)-1KH=KU-KLIF(KH)304,240,210210 K=N-KHIC=0KLT=KUDO 260 J=1,KHKLT=KLT-1IC=IC+1KI=MAXA(K)ND=MAXA(K+1)-KI-1IF(ND) 260,260,270270 KK=MIN0(IC,ND)C=0.0DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L)A(KLT)=A(KLT)-C260 K=K+1240 K=NB=0.0DO 300 KK=KL,KUK=K-1KI=MAXA(K)C=A(KK)/A(KI)IF(ABS(C).LT.1.0E+07) GOTO 290WRITE(IOUT,2010) N,CSTOP290 B=B+C*A(KK)300 A(KK)=CA(KN)=A(KN)-B304 IF(A(KN)) 310,310,200310 IF(ISH.EQ.0) GOTO 320IF(A(KN).EQ.0.0) A(KN)=-1.0E-16GOTO 200320 WRITE(IOUT,2000) N,A(KN)STOP200 CONTINUERETURN2000 FORMAT(//,' Stop-stiffness matrix not positive definite',//,'no positive& pivot for equation',I4,&//,' pivot =',E20.10)2010 FORMAT(//,' Stop-strum sequence check failed + because of multiplier& growth for column &number',I4,//,' Multiplier = ',E20.8)END!回代，求得节点位移SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM)DIMENSION A(NWK),V(NN,1),MAXA(NNM)NIP=1DO IP=1,NIPDO 400 N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 400,410,410410 K=NC=0.0DO 420 KK=KL,KUK=K-1420 C=C+A(KK)*V(K,IP)V(N,IP)=V(N,IP)-C400 CONTINUEDO 480 N=1,NNK=MAXA(N)480 V(N,IP)=V(N,IP)/A(K)IF(NN.EQ.1)RETURNN=NNDO 500 L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 500,510,510510 K=NDO 520 KK=KL,KUK=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1ENDDORETURNEND!求解杆件内力、支反力和位移SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y, Z,PP,FF,SG,SM,ME,& NR,RR,NF)DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),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),PP(NPF),&SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)SG=0;SM=0;FF=0;FF2=0DO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0) THENDIST(3*(I-1)+J)=0.0ELSEIF(B.LE.NN) THENDIST(3*(I-1)+J)=FTOOL(LAB)ENDIFENDDOENDDODO IE=1,NEN1=ME(1,IE);N2=ME(2,IE)DO J=1,NFUE(J)=DIST(3*(N1-1)+J)UE(3+J)=DIST(3*(N2-1)+J)ENDDOCALL 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,NFFF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOISW=NAE(IE)AO=AE(2,ISW)SG(IE)=FE1(2)SM(IE)=FE1(2)/AODO I=1,NPFFF2(I)=FF(I)-PP(I)ENDDOENDDODO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)' 'WRITE(12,*)'****************************************'WRITE(12,*)'*********The Results of Calculation**********'WRITE(12,*)'****************************************'WRITE(12,600)WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000, I=1,NP)WRITE(12,620)WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)WRITE(12,640)WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)600 FORMAT(6X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)') 610 FORMAT(1X,I4,2X,1P3E12.2)620 FORMAT(//6X,'The Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'σ(MPa)')630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)640 FORMAT(//6X,'The Bearing Force'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')650 FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（b）所示。
桁架结构MATLAB编程
桁架结构% 变量说明% NPOIN NELEM NVFIX NFPOIN NFPRES% 总结点数,单元数, 约束个数, 受力结点数, 非结点力数% COORD LNODS YOUNG% 结构节点坐标数组, 单元定义数组, 弹性模量% FPOIN FPRES FORCE FIXED% 结点力数组，非结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵,结点位移向量%**************************************************format short e%设定输出类型clear %清除内存变量FP1=fopen('6-6.txt','rt'); %打开初始数据文件%读入控制数据NELEM=fscanf(FP1,'%d',1); %单元数NPOIN=fscanf(FP1,'%d',1); %结点数NVFIX=fscanf(FP1,'%d',1); %约束数NFPOIN=fscanf(FP1,'%d',1); %结点荷载数NFPRES=fscanf(FP1,'%d',1); %非结点荷载数YOUNG=fscanf(FP1,'%f',1); %弹性模量% 读取结构信息LNODS=fscanf(FP1,'%f',[4,NELEM])'% 单元定义：左、右结点号，面积，惯性矩（共计 NELEM组）COORD=fscanf(FP1,'%f',[2,NPOIN])'% 坐标： x,y坐标（共计 NPOIN 组）FPOIN=fscanf(FP1,'%f',[3,NFPOIN])'% 节点力（共计 NFPOIN 组）：结点号、X方向力（向右正），% Y方向力（向上正），FPRES=fscanf(FP1,'%f',[4,NFPRES])' % 均布力（共计% NFPRES 组）：单元号、荷载类型、荷载大小、距离左端长度FIXED=fscanf(FP1,'%f',NVFIX)'% 约束信息：约束对应的位移编码（共计 NVFIX 组）%---------------------------------------------------------HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零%形成总刚for i=1:NELEM % 对单元个数循环% 生成局部单刚(局部坐标) 右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD); % 坐标转换矩阵EKT=T'*EK*T; % 生成整体单刚(整体坐标系)% 组成总刚按2*2子块加入总刚中（共计4块)for j=1:2 %对行进行循环---按结点号循环N1=LNODS(i,j)*2; % j结点第2个位移的整体编码for k=1:2 %对列进行循环---按结点号循环N2=LNODS(i,k)*2; % k结点第2个位移的整体编码HK((N1-1):N1,(N2-1):N2)=HK((N1-1):N1,(N2-1):N2)...+EKT(j*2-1:j*2,k*2-1:k*2);endendend% 由结点力与非结点力生成总荷载向量列阵for i=1:NFPOIN % 对结点荷载个数进行循环N1=FPOIN(i,1); % 作用荷载的结点号N1=N1*2-2; % 该结点号对应第一个位移编码 - 1for j=1:2FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取结点荷载endend% 计算由非结点荷载引起的等效结点荷载for i=1:NFPRES % 对非结点荷载个数进行循环F0=ele_FPRES(i,FPRES,LNODS,COORD); %计算单元固端力% 对单元局部杆端力要进行坐标转换T=zbzh(i,LNODS,COORD); % 坐标转换矩阵F0=T'*F0;ele=FPRES(i,1); % 取荷载所在的单元号NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号% 将单元固端力变成等效结点荷载（注意固端力与等效结点荷载符号相反）FORCE((3*NL-2):3*NL)=FORCE((3*NL-2):3*NL)-F0(1:3);FORCE((3*NR-2):3*NR)=FORCE((3*NR-2):3*NR)-F0(4:6);end% 总刚、总荷载进行边界条件处理for j=1:NVFIX % 对约束个数进行循环N1=FIXED(j);HK(1:2*NPOIN,N1)=0; HK(N1,1:2*NPOIN)=0; HK(N1,N1)=1;% 将零位移约束对应的行、列变成零，主元变成1FORCE(N1)=0;end%---------------------------------------------------------DISP=HK\FORCE % 方程求解，HK先求逆再与力向量左乘%---------------------------------------------------------% 求结构各个单元内力EDISP=zeros(4,1); % 单元位移列向量清零for i=1:NELEM % 对单元个数进行循环for j=1:2 %对杆端循环% i单元左右端结点号*2 = 该结点的最后一个位移编码N1=LNODS(i,j)*2;% 取一端的单元位移列向量EDISP(2*j-1:2*j)=DISP(N1-1:N1);end% 生成局部单刚(局部坐标) 右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD); % 坐标转换矩阵FE=EK*T*EDISP; %计算局部坐标杆端力（由结点位移产生）for j=1:NFPRESif FPRES(j,1) == i %成立时，当前单元上有非结点荷载F0=ele_FPRES(j,FPRES,LNODS,COORD);%单元固端力FE=FE+F0; % 考虑由非结点荷载引起的杆端力endendFE % 打印杆端力end%-------------------------------------------------------------------------------- % 计算桁架单元刚度矩阵函数 EK% 入口参数：单元号、单元信息数组、结点坐标、弹性模量% 出口参数：局部桁架单元刚度矩阵EKfunction EK=ele_EK(i,LNODS,COORD,YOUNG)NL=LNODS(i,1); NR=LNODS(i,2); %左右结点号dx=COORD(NR,1)-COORD(NL,1); % x 坐标差dy=COORD(NR,2)-COORD(NL,2); % y 坐标差L=sqrt(dx^2+dy^2); %单元长度A=LNODS(i,3); %面积；E=YOUNG;% 生成单刚(局部坐标) 右手坐标系EK =[E*A/L 0 -E*A/L 0 ;...0 0 0 0 ;...-E*A/L 0 E*A/L 0 ;...0 0 0 0 ];return%---------------------------------------------------------------------------------%计算单元固端力函数（正方向：X向右 Y向上 M逆时针）% 入口参数：荷载序号，荷载信息，单元信息，结点坐标% 出口参数：单元固端力——左右两端的轴力、剪力、弯矩function F0=ele_FPRES(iFPRES,FPRES,LNODS,COORD)ele=FPRES(iFPRES,1); %取荷载所在的单元号G=FPRES(iFPRES,3); %单元荷载大小C=FPRES(iFPRES,4); %单元荷载与左端距离NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号dx=COORD(NR,1)-COORD(NL,1); % x 坐标差dy=COORD(NR,2)-COORD(NL,2); % y 坐标差L=sqrt(dx^2+dy^2); %单元长度% 计算公式中一些常出现的项D=L-C; C1=C/L; C2=C1*C1; C3=C1*C2;B1=D/L; B2=B1/L;F0=[0;0;0;0;0;0]; %单元固端力清零switch FPRES(iFPRES,2)case 1 %均布荷载F0(2)=-G*C*(2-2*C2+C3)/2.0;F0(3)=-G*C*C*(6-8*C1+3*C2)/12.0;F0(5)=-G*C-F0(2);F0(6)=G*C*C*C1*(4-3*C1)/12.0;case 2 %横向集中力F0(2)=-G*B1*B2*(L+2*C);F0(3)=-G*C*B1*B1;F0(5)=-G*C2*(L+2*D)/L;F0(6)=G*D*C2;case 3 %纵向集中力F0(1)=-G*B1;F0(4)=-G*C1;endreturn%-------------------------------------------------------------------------------- % 形成第i单元的坐标转换矩阵函数 T(4,4)% 入口参数：单元号，单元信息，结点坐标% 出口参数：坐标转换矩阵（整体向局部投影）function T=zbzh(i,LNODS,COORD)NL=LNODS(i,1); %左结点号NR=LNODS(i,2); %右结点号dx=COORD(NR,1)-COORD(NL,1); % x 坐标差dy=COORD(NR,2)-COORD(NL,2); % y 坐标差L=sqrt(dx^2+dy^2); % 单元长度c=dx/L; % cos a (与 x 轴夹角余弦） s=dy/L; % sin aT=[ c s 0 0;...-s c 0 0;...0 0 c s0 0 -s c ];return6-6.txt文件数据：4 45 2 0 2.95e111 2 0.001 13 2 0.001 11 3 0.001 14 3 0.001 10 00.4 00.4 0.30 0.32 2e4 03 0 -2.5e41 2 4 7 8。
新代桁架系统编程教程
新代桁架系统编程教程引言：新代桁架系统是一种新兴的结构系统，通过桁架结构的优势，能够在建筑领域中发挥重要作用。
而对于新代桁架系统的编程，更是为设计师和工程师们提供了更多的可能性。
本教程将以新代桁架系统编程为主题，为读者详细介绍如何进行新代桁架系统的编程。
一、新代桁架系统的概述新代桁架系统是一种基于传统桁架结构的升级改进，它采用了更加先进的设计理念和工艺技术。
新代桁架系统在结构强度、搭建速度和灵活性方面都有显著的提升，因此在建筑领域中得到了广泛应用。
而新代桁架系统的编程则是为了更好地实现其功能和优势。
二、新代桁架系统编程的基础知识1. 桁架结构的原理和特点：了解桁架结构的基本原理和特点，包括节点、杆件、受力分析等内容，为后续的编程打下基础。
2. 编程软件的选择：选择合适的编程软件，例如Rhino、Grasshopper等，以及相应的插件或脚本工具，为编程工作提供支持。
3. 编程语言的学习：学习一种适合新代桁架系统编程的语言，例如Python、C#等，掌握其基本语法和常用函数，为编程工作做好准备。
三、新代桁架系统编程的实践步骤1. 数据准备：收集和整理与新代桁架系统相关的数据，包括结构参数、材料属性、受力要求等，为后续的编程工作提供输入。
2. 桁架生成：通过编程软件，根据输入的数据和设计要求，生成新代桁架系统的初始结构模型，并进行初步的受力分析。
3. 结构优化：通过编程工具，对新代桁架系统进行优化设计，以提高结构的性能和效率，例如减少材料使用量、优化节点连接等。
4. 受力分析：利用编程软件进行新代桁架系统的受力分析，验证结构的稳定性和强度，确保其满足设计要求。
5. 结果展示：通过编程软件将优化后的新代桁架系统进行可视化展示，以便设计师和工程师们更好地理解和评估结构的性能。
四、新代桁架系统编程的应用案例1. 建筑结构设计：通过新代桁架系统的编程，实现建筑结构的设计和优化，提高结构的稳定性和安全性。
2. 建筑外观设计：利用新代桁架系统的编程，实现建筑外观的设计和变形效果，为建筑增添独特的艺术魅力。
空间桁架结构程序的设计(Fortran)
空间桁架静力分析程序及算例1、变量及数组说明2、空间桁架结构有限元分析程序源代码!主程序（读入文件，调用总计算程序，输出结果）CHARACTER IDFUT*20,OUTFUT*20WRITE(*,*) 'Input Data File name:'READ (*,*)IDFUTOPEN (11,FILE=IDFUT,STATUS='OLD')WRITE(*,*) 'Output File name:'READ (*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')WRITE(12,*)'*****************************************'WRITE(12,*)'* Program for Analysis of Space Trusses *'WRITE(12,*)'* School of Civil Engineering CSU *'WRITE(12,*)'* 2012.6.25 Designed By MuZhaoxiang *'WRITE(12,*)'*****************************************'WRITE(12,*)' 'WRITE(12,*)'*****************************************'WRITE(12,*)'*************The Input Data****************'WRITE(12,*)'*****************************************'WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100 FORMAT(6X,'The General Information'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',& 5X,'NCF',5X,'ND')110 FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!********************************************************************!总计算程序SUBROUTINE ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR), NAE(NE),&AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),&PP(NPF),FF(NPF),SG(NE),SM(NE)READ(11,*)(X(I),Y(I),Z(I),I=1,NP)READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)READ(11,*)(RR(1,J),RR(2,J),J=1,NR)READ(11,*)(AE(1,J),J=1,2)WRITE(12,120)WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)WRITE(12,130)WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)WRITE(12,140)WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)WRITE(12,150)WRITE(12,151)(AE(1,J),J=1,2)IF(NCF/=0)THENREAD(11,*)((PF(I,J),I=1,4),J=1,NCF)WRITE(12,160)WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)ENDIF120 FORMAT(/6X,'The Information of Joints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')121 FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,'The Information of Members'/2x,'Member',2X,'START',4X,'END',6X,'NAE')131 FORMAT(1X,I4,3I8)140 FORMAT(/6X,'The Information of SUPPORTS'/2x,'Joint',5X,'S')141 FORMAT(1X,I4,F8.3)150 FORMAT(/6X,'The Information of Sections'/4x,'E0',8X,'A0')151 FORMAT(1X,1PE8.2,F8.4)160 FORMAT(/6X,'The Loading at Joints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')161 FORMAT(1X,I4,3F8.2)CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)ISH=1CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALL REBACK(CKK,V,MAXA,NN,NWK,NNM)CALL DISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF)END!********************************************************************!矩阵转置子程序SUBROUTINE MAT(M,N,A,B)DIMENSION A(M,N),B(N,M)DO I=1,MDO J=1,NB(J,I)=A(I,J)END DOEND DORETURNEND!单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)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/BLAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND!单元坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)T=0N1=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)/BLCY=(Y2-Y1)/BLCZ=(Z2-Z1)/BLT(1,1)=CX;T(2,4)=CXT(1,2)=CY;T(2,5)=CYT(1,3)=CZ;T(2,6)=CZRETURNEND!生成单元联系数组LMTSUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0DO I=1,NPC=0DO K=1,NRKR=RR(1,K)IF(KR.EQ.I) C=RR(2,K)ENDDONC=C !NC=0,提取了整数部分C=C-NC !C=0.***,例如C=0.111DO J=1,NFC=C*10.0 !例如C=1.21L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位 !上的数字,这里"+0.1"是为了防止四舍五入是出现错误 C=C-LIF(L.EQ.0)THENN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DO IE=1,NEDO I=1,NDNI=ME(I,IE)DO J=1,NFLMT((I-1)*NF+J,IE)=IT(J,NI)ENDDOENDDOENDDORETURNEND!二维总刚中对角线元地址数组SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)DIMENSION MAXA(NPF),LMT(NDF,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NNMIP=I-1IG=IPDO IE=1,NEDO J=1,NDFIF(LMT(J,IE).EQ.IP) THENDO K=1,NDFIF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE)ENDDOEND IFENDDOENDDOMAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK= MAXA(NNM)-1RETURNEND!生成一维存储结构总刚度矩阵SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),& MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)CKK=0DO 10 IE=1,NETAK=0CALL 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,6DO 220 J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IF((NJ-NI).GE.0.AND.NI*NJ.GT.0) THENIJ=MAXA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIF220 CONTINUE10 CONTINUERETURNEND!生成荷载矩阵SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF)V=0PP=0DO I=1,NFDO J=1,NPDO K=1,NCFIF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J))=PF(I+1,K)ENDIFENDDOENDDOENDDODO K=1,NCFDO I=1,NPIF(I.EQ.PF(1,K))THENPP(NF*(I-1)+1)=PF(2,K)PP(NF*(I-1)+2)=PF(3,K)PP(NF*(I-1)+3)=PF(4,K)ENDIFENDDOENDDORETURNEND!对一维结构总刚度矩阵进行矩阵分解（LDLT）SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM) DIMENSION A(NWK),MAXA(NNM)IF(NN.EQ.1) RETURNDO 200 N=1,NNKN=MAXA(N)KL=KN+1KU=MAXA(N+1)-1KH=KU-KLIF(KH)304,240,210210 K=N-KHIC=0KLT=KUDO 260 J=1,KHKLT=KLT-1IC=IC+1KI=MAXA(K)ND=MAXA(K+1)-KI-1IF(ND) 260,260,270270 KK=MIN0(IC,ND)C=0.0DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L)A(KLT)=A(KLT)-C260 K=K+1240 K=NB=0.0DO 300 KK=KL,KUK=K-1KI=MAXA(K)C=A(KK)/A(KI)IF(ABS(C).LT.1.0E+07) GOTO 290WRITE(IOUT,2010) N,CSTOP290 B=B+C*A(KK)300 A(KK)=CA(KN)=A(KN)-B304 IF(A(KN)) 310,310,200310 IF(ISH.EQ.0) GOTO 320IF(A(KN).EQ.0.0) A(KN)=-1.0E-16GOTO 200320 WRITE(IOUT,2000) N,A(KN)STOP200 CONTINUERETURN2000 FORMAT(//,' Stop-stiffness matrix not positive definite',//,'no positive& pivot for equation',I4,&//,' pivot =',E20.10)2010 FORMAT(//,' Stop-strum sequence check failed + because of multiplier& growth for column &number',I4,//,' Multiplier = ',E20.8)END!回代，求得节点位移SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM)DIMENSION A(NWK),V(NN,1),MAXA(NNM)NIP=1DO IP=1,NIPDO 400 N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 400,410,410410 K=NC=0.0DO 420 KK=KL,KUK=K-1420 C=C+A(KK)*V(K,IP)V(N,IP)=V(N,IP)-C400 CONTINUEDO 480 N=1,NNK=MAXA(N)480 V(N,IP)=V(N,IP)/A(K)IF(NN.EQ.1)RETURNN=NNDO 500 L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 500,510,510510 K=NDO 520 KK=KL,KUK=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1ENDDORETURNEND!求解杆件力、支反力和位移SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y, Z,PP,FF,SG,SM,ME,& NR,RR,NF)DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),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),PP(NPF),&SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)SG=0;SM=0;FF=0;FF2=0DO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0) THENDIST(3*(I-1)+J)=0.0ELSEIF(B.LE.NN) THENDIST(3*(I-1)+J)=FTOOL(LAB)ENDIFENDDOENDDODO IE=1,NEN1=ME(1,IE);N2=ME(2,IE)DO J=1,NFUE(J)=DIST(3*(N1-1)+J)UE(3+J)=DIST(3*(N2-1)+J)ENDDOCALL 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,NFFF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOISW=NAE(IE)AO=AE(2,ISW)SG(IE)=FE1(2)SM(IE)=FE1(2)/AODO I=1,NPFFF2(I)=FF(I)-PP(I)ENDDOENDDODO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)' 'WRITE(12,*)'****************************************'WRITE(12,*)'*********The Results of Calculation**********'WRITE(12,*)'****************************************'WRITE(12,600)WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000, I=1,NP)WRITE(12,620)WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)WRITE(12,640)WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)600 FORMAT(6X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)') 610 FORMAT(1X,I4,2X,1P3E12.2)620 FORMAT(//6X,'T he Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'σ(MPa)')630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)640 FORMAT(//6X,'The Bearing Force'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')650 FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（b）所示。
空间桁架程序设计
第六章 空间桁架程序设计第一节 概述一、计算模型及计算方法1． 计算模型及整体坐标系下单元刚度矩阵的形成 选取计算模型时，应以杆简联结点和支座结点作为计算结点，各结点均为光滑的理想铰结点；以任意两结点间的杆简为计算单元，各单元支承受轴力；非结点荷载要转化为等效结点荷载，各单元抗拉压刚度相同。
由于结构为空间桁架，所以，结构各结点的位移分量为[]Tz y x ∆∆∆=∆在局部坐标系下，单元的杆端位移列阵和杆端力列阵分别为e j j j i i i ej i ew v u w v u ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎦⎤⎢⎢⎣⎡=δδδ ej j j i i i e j i e z y x z y x F F F ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎦⎤⎢⎢⎣⎡= （6-1） 式中j j j i i i w v u w v u 、、和、、分别为结点i 、j 沿局部坐标系z y x 、、方向的线位移，见图6-1。
杆端力应与杆端位移一一对应，图中没再绘出。
当单元的杆端位移分量为任意值时，可写出空间桁架单元刚度方程。
以矩阵表示为ej j j i i i eej j j i i i w v u w v u l EA Z Y X Z Y X ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡00000000001001000000000000001001 （6-2）简写成 ee e k F δ= （6-3） 式中eel EA k ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=000000000001001000000000000001001 （6-4） 称为单元○e 在局部坐标系下的刚度矩阵。
在空间桁架中，各杆方向不尽相同。
每根杆件采用各自的局部坐标系，这对于单元分析十分方便。
空间桁架结构程序设计(Fortran)
空间桁架静力分析程序及算例1、变量及数组说明2、空间桁架结构有限元分析程序源代码!主程序（读入文件，调用总计算程序，输出结果）CHARACTER IDFUT*20,OUTFUT*20WRITE(*,*) 'Input Data File name:'READ (*,*)IDFUTOPEN (11,FILE=IDFUT,STATUS='OLD')WRITE(*,*) 'Output File name:'READ (*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')WRITE(12,*)'*****************************************'WRITE(12,*)'* Program for Analysis of Space Trusses *'WRITE(12,*)'* School of Civil Engineering CSU *'WRITE(12,*)'* 2012.6.25 Designed By MuZhaoxiang *'WRITE(12,*)'*****************************************'WRITE(12,*)' 'WRITE(12,*)'*****************************************'WRITE(12,*)'*************The Input Data****************'WRITE(12,*)'*****************************************'WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100 FORMAT(6X,'The General Information'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',& 5X,'NCF',5X,'ND')110 FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!********************************************************************!总计算程序DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR), NAE(NE),& AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),&PP(NPF),FF(NPF),SG(NE),SM(NE)READ(11,*)(X(I),Y(I),Z(I),I=1,NP)READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)READ(11,*)(RR(1,J),RR(2,J),J=1,NR)READ(11,*)(AE(1,J),J=1,2)WRITE(12,120)WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)WRITE(12,130)WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)WRITE(12,140)WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)WRITE(12,150)WRITE(12,151)(AE(1,J),J=1,2)IF(NCF/=0)THENREAD(11,*)((PF(I,J),I=1,4),J=1,NCF)WRITE(12,160)WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)ENDIF120 FORMAT(/6X,'The Information of Joints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')121 FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,'The Information of Members'/2x,'Member',2X,'START',4X,'END',6X,'NAE')131 FORMAT(1X,I4,3I8)140 FORMAT(/6X,'The Information of SUPPORTS'/2x,'Joint',5X,'S')141 FORMAT(1X,I4,F8.3)150 FORMAT(/6X,'The Information of Sections'/4x,'E0',8X,'A0')151 FORMAT(1X,1PE8.2,F8.4)160 FORMAT(/6X,'The Loading at Joints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')161 FORMAT(1X,I4,3F8.2)CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)ISH=1CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALL REBACK(CKK,V,MAXA,NN,NWK,NNM)CALL DISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF)END!********************************************************************!矩阵转置子程序SUBROUTINE MAT(M,N,A,B)DIMENSION A(M,N),B(N,M)DO I=1,MDO J=1,NEND DOEND DORETURNEND!单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)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/BLAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND!单元坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)T=0N1=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)/BLCY=(Y2-Y1)/BLCZ=(Z2-Z1)/BLT(1,1)=CX;T(2,4)=CXT(1,2)=CY;T(2,5)=CYT(1,3)=CZ;T(2,6)=CZRETURNEND!生成单元联系数组LMTSUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT) DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0C=0DO K=1,NRKR=RR(1,K)IF(KR.EQ.I) C=RR(2,K)ENDDONC=C !NC=0,提取了整数部分C=C-NC !C=0.***,例如C=0.111DO J=1,NFC=C*10.0 !例如C=1.21L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位 !上的数字,这里"+0.1"是为了防止四舍五入是出现错误C=C-LIF(L.EQ.0)THENN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DO IE=1,NEDO I=1,NDNI=ME(I,IE)DO J=1,NFLMT((I-1)*NF+J,IE)=IT(J,NI)ENDDOENDDOENDDORETURNEND!二维总刚中对角线元地址数组SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)DIMENSION MAXA(NPF),LMT(NDF,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NNMIP=I-1IG=IPDO IE=1,NEDO J=1,NDFIF(LMT(J,IE).EQ.IP) THENDO K=1,NDFENDDOEND IFENDDOENDDOMAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK= MAXA(NNM)-1RETURNEND!生成一维存储结构总刚度矩阵SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM) DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),& MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)CKK=0DO 10 IE=1,NETAK=0CALL 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,6DO 220 J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IF((NJ-NI).GE.0.AND.NI*NJ.GT.0) THENIJ=MAXA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIF220 CONTINUE10 CONTINUERETURNEND!生成荷载矩阵SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF)V=0PP=0DO I=1,NFDO J=1,NPDO K=1,NCFIF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J))=PF(I+1,K)ENDIFENDDOENDDODO K=1,NCFDO I=1,NPIF(I.EQ.PF(1,K))THENPP(NF*(I-1)+1)=PF(2,K)PP(NF*(I-1)+2)=PF(3,K)PP(NF*(I-1)+3)=PF(4,K)ENDIFENDDOENDDORETURNEND!对一维结构总刚度矩阵进行矩阵分解（LDLT）SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM)DIMENSION A(NWK),MAXA(NNM)IF(NN.EQ.1) RETURNDO 200 N=1,NNKN=MAXA(N)KL=KN+1KU=MAXA(N+1)-1KH=KU-KLIF(KH)304,240,210210 K=N-KHIC=0KLT=KUDO 260 J=1,KHKLT=KLT-1IC=IC+1KI=MAXA(K)ND=MAXA(K+1)-KI-1IF(ND) 260,260,270270 KK=MIN0(IC,ND)C=0.0DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L)A(KLT)=A(KLT)-C260 K=K+1240 K=NB=0.0DO 300 KK=KL,KUK=K-1KI=MAXA(K)C=A(KK)/A(KI)WRITE(IOUT,2010) N,CSTOP290 B=B+C*A(KK)300 A(KK)=CA(KN)=A(KN)-B304 IF(A(KN)) 310,310,200310 IF(ISH.EQ.0) GOTO 320IF(A(KN).EQ.0.0) A(KN)=-1.0E-16GOTO 200320 WRITE(IOUT,2000) N,A(KN)STOP200 CONTINUERETURN2000 FORMAT(//,' Stop-stiffness matrix not positive definite',//,'no positive& pivot for equation',I4,& //,' pivot =',E20.10)2010 FORMAT(//,' Stop-strum sequence check failed + because of multiplier& growth for column & number',I4,//,' Multiplier = ',E20.8)END!回代，求得节点位移SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM)DIMENSION A(NWK),V(NN,1),MAXA(NNM)NIP=1DO IP=1,NIPDO 400 N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 400,410,410410 K=NC=0.0DO 420 KK=KL,KUK=K-1420 C=C+A(KK)*V(K,IP)V(N,IP)=V(N,IP)-C400 CONTINUEDO 480 N=1,NNK=MAXA(N)480 V(N,IP)=V(N,IP)/A(K)IF(NN.EQ.1)RETURNN=NNDO 500 L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 500,510,510510 K=NDO 520 KK=KL,KU520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1ENDDORETURNEND!求解杆件内力、支反力和位移SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y, Z,PP,FF,SG,SM,ME,& NR,RR,NF)DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),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),PP(NPF),&SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)SG=0;SM=0;FF=0;FF2=0DO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0) THENDIST(3*(I-1)+J)=0.0ELSEIF(B.LE.NN) THENDIST(3*(I-1)+J)=FTOOL(LAB)ENDIFENDDOENDDODO IE=1,NEN1=ME(1,IE);N2=ME(2,IE)DO J=1,NFUE(J)=DIST(3*(N1-1)+J)UE(3+J)=DIST(3*(N2-1)+J)ENDDOCALL 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,NFFF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOISW=NAE(IE)AO=AE(2,ISW)SG(IE)=FE1(2)SM(IE)=FE1(2)/AODO I=1,NPFFF2(I)=FF(I)-PP(I)ENDDODO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)' 'WRITE(12,*)'****************************************'WRITE(12,*)'*********The Results of Calculation**********'WRITE(12,*)'****************************************'WRITE(12,600)WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000, I=1,NP)WRITE(12,620)WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)WRITE(12,640)WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)600 FORMAT(6X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)')610 FORMAT(1X,I4,2X,1P3E12.2)620 FORMAT(//6X,'The Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'σ(MPa)')630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)640 FORMAT(//6X,'The Bearing Force'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')650 FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（b）所示。
3d3s空间桁架设计步骤 理论说明以及概述
3d3s空间桁架设计步骤理论说明以及概述1. 引言1.1 概述本文将详细介绍3D3S空间桁架设计的步骤、理论说明以及概述。
空间桁架是一种由杆件和节点组成的结构体系，具有轻质、高强度和刚性好等特点，并被广泛应用于建筑工程、航空航天领域以及体育场馆等各个领域。
在本文中，我们将对3D3S空间桁架设计的过程进行深入研究，并探讨其相关的理论知识。
1.2 文章结构本文共分为四个部分。
首先，在引言部分，我们将给出文章的概述，明确文章结构以及研究目的。
然后，在第二部分中，我们将详细介绍3D3S空间桁架设计的步骤，包括确定设计需求、建立初始几何模型以及分析和优化设计。
接下来，在第三部分中，我们将对空间桁架原理进行概述，并介绍与之相关的结构力学基础知识以及桁架结构的特点和应用领域。
最后，在结论部分，我们将总结本文所介绍的设计步骤和理论知识，并展望未来的发展方向。
1.3 目的本文的目的是为读者提供关于3D3S空间桁架设计步骤、理论说明以及概述的全面了解。
通过本文的阐述，读者将能够掌握从确定设计需求到最终优化设计的全过程，并能够理解空间桁架原理以及相关的结构力学基础知识。
同时，本文也旨在激发读者对未来空间桁架设计发展方向的思考，并为相关领域的专业人士提供参考和借鉴。
通过深入研究和探索，我们相信这篇文章将对3D3S空间桁架设计领域有所贡献。
2. 3D3S空间桁架设计步骤:2.1 确定设计需求:在进行3D3S（三维数字结构的静力分析与设计）空间桁架设计之前，需要首先明确设计的具体需求。
这包括结构用途、预期负荷、支撑和连接要求等方面的考虑。
通过确定设计需求，可以为后续的步骤提供清晰的目标和指导。
2.2 建立初始几何模型:在开始详细设计之前，需要建立一个初始的几何模型。
这一步骤包括确定桁架的整体形状、尺寸和布局。
可以使用专业CAD软件或者手工绘图来创建这个初步模型。
2.3 分析和优化设计:接下来是对初始几何模型进行分析和优化。
桁架结构 练习5收敛准则
5.收敛准则 …桁架结构
6.加紧力收敛准则 Main Menu > Solution > Solve > Sol’n Control 选择表：“Nonlinear” 设置收敛准则 … 替换 … 结构力 F TOLER 容差大概值 = 0.0001 NORM 收敛范数 = L2 norm MINREF 最小参考值 = 1.0 [确定] [关闭] [确定] 或输入如下命令： /SOLUTION CNVTOL, F, , 0.0001, , 1.0
4
4.图示垂直位移 Main Menu > General Postproc > Plot Results > Contour Plot > Nodal Solu 自由度解 Y 向位移 [确定] 或输入如下命令： /POST1 PLNSOL, U, Y 注意此分析的位移最大值： DMX = 4.1 inches 5.查看监控文件的内容 (truss.mntr)
数据点 # 1
数据点 # 2
数据点 # 3
输入的应力-应变曲线
数据点 # 1
数据点 # 2
缺省准则
迭代#2
迭代 #1
牛顿-拉斐逊迭代
施加的 应力
收敛解
5.收敛准则 …桁架结构
5.收敛准则 …桁架结构
5.收敛准则 …桁架结构
解决此非线性塑性分析（用缺省的收敛控制）
Main Menu > Solution > Solve > Current LS [确定]
2
或输入如下命令： /SOLUTION /STATUS, SOLUTION SOLVE
3
注意：当前收敛控制设为缺省设置。
注: 当前收敛控制设置为用户指定值.
桁架资料_精品文档
多边形桁架也称折线形桁架(图c)。
上弦节点位于二次抛物线上，如上弦呈拱形可减少节间荷载产生的弯矩,但制造较为复杂。
在均布荷载作用下,桁架外形和简支梁的弯矩图形相似，于是上下弦轴力分布均匀，腹杆轴力较小，用料最省，是工程中常用的一种桁架形式。
形式选择 从力学方面分析，桁架外形与简支梁的弯矩图相似时，上下弦杆的轴力分布均匀，腹杆轴力小，用料最省；从材料与创造方面分析，木桁架做成三角形，钢桁架采用梯形或者平行弦形，钢筋混凝土与预应力混凝土桁架为多边形或者梯形为宜。
平面桁架普通按理想的铰接桁架进行计算，即假设荷载施加在桁架节点上（如果荷载施加在节间时，可按简支梁换算为节点荷载），并和桁架的全部杆件均在同一平面内，杆件的重心轴在向来线上，节点为可自由转动的铰接点。
理想状态下的静定桁架，可以将杆件轴力作为未知量，按静力学的数解法或者图解法求出已知荷载下杆件的轴向拉力或者压力工程用的桁架节点，普通是具有一定刚性的节点而不是理想的铰接节点，由于节点刚性的影响而浮现的杆件弯曲应力和轴向应力称为次应力。
计算次应力需考虑杆件轴向变形，可用超静定结构的方法或者有限元法求解。
空间桁架由若干个平面桁架所组成，可将荷载分解成与桁架同一平面的分力按平面桁架进行计算，或者按空间铰接杆系用有限元法计算木桁架设计原则桁架的高跨比1.桁架跨度中央的高度h与跨度l的比值称为高跨比。
为保证桁架具有足够的刚度，按桁架的外形，分别规所示。
定木桁架、钢木桁架高跨比的最小限值如表6.3.1规定的桁架，不必再核算其挠度。
高跨比值已符合表6.3.1桁架的预起拱度2.为了消除桁架可见的挠度，不论木桁架或者钢木桁架，皆应在创造时预先向上起拱。
起拱度通常取为桁架跨。
起拱时应保持桁架的高跨比不变，木桁架常在下弦接头处提高(图6.3.9)，而钢木桁架则常在度的1/200下弦节点处提高。
桁架节间的划分3.桁架节间的划分原则是：根据荷载、跨度及所用木材强度设计值的大小进行节间划分，在常用木材规格范围内，充分利用上弦的承载能力。
空间桁架结构程序设计(Fortran)
空间桁架静力分析程序及算例1、变量及数组说明2、空间桁架结构有限元分析程序源代码!主程序（读入文件，调用总计算程序，输出结果）CHARACTER IDFUT*20,OUTFUT*20WRITE(*,*) 'Input Data File name:'READ (*,*)IDFUTOPEN (11,FILE=IDFUT,STATUS='OLD')WRITE(*,*) 'Output File name:'READ (*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')WRITE(12,*)'*****************************************'WRITE(12,*)'* Program for Analysis of Space Trusses *'WRITE(12,*)'* School of Civil Engineering CSU *'WRITE(12,*)'* 2012.6.25 Designed By MuZhaoxiang *'WRITE(12,*)'*****************************************'WRITE(12,*)' 'WRITE(12,*)'*****************************************'WRITE(12,*)'*************The Input Data****************'WRITE(12,*)'*****************************************'WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100 FORMAT(6X,'The General Information'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',&5X,'NCF',5X,'ND')110 FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!******************************************************************** !总计算程序SUBROUTINE ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR), NAE(NE),&AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),&PP(NPF),FF(NPF),SG(NE),SM(NE)READ(11,*)(X(I),Y(I),Z(I),I=1,NP)READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)READ(11,*)(RR(1,J),RR(2,J),J=1,NR)READ(11,*)(AE(1,J),J=1,2)WRITE(12,120)WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)WRITE(12,130)WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)WRITE(12,140)WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)WRITE(12,150)WRITE(12,151)(AE(1,J),J=1,2)IF(NCF/=0)THENREAD(11,*)((PF(I,J),I=1,4),J=1,NCF)WRITE(12,160)WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)ENDIF120 FORMAT(/6X,'The Information of Joints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z') 121 FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,'The Information of Members'/2x,'Member',2X,'START',4X,'END',6X,'NAE')131 FORMAT(1X,I4,3I8)140 FORMAT(/6X,'The Information of SUPPORTS'/2x,'Joint',5X,'S')141 FORMAT(1X,I4,F8.3)150 FORMAT(/6X,'The Information of Sections'/4x,'E0',8X,'A0')151 FORMAT(1X,1PE8.2,F8.4)160 FORMAT(/6X,'The Loading at Joints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ') 161 FORMAT(1X,I4,3F8.2)CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)ISH=1CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALL REBACK(CKK,V,MAXA,NN,NWK,NNM)CALL DISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF) END!********************************************************************!矩阵转置子程序SUBROUTINE MAT(M,N,A,B)DIMENSION A(M,N),B(N,M) DO I=1,MDO J=1,NB(J,I)=A(I,J)END DOEND DORETURNEND!单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)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/BLAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND!单元坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)T=0N1=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)/BLCY=(Y2-Y1)/BLCZ=(Z2-Z1)/BLT(1,1)=CX;T(2,4)=CXT(1,2)=CY;T(2,5)=CYT(1,3)=CZ;T(2,6)=CZRETURNEND!生成单元了解数组LMTSUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0DO I=1,NPC=0DO K=1,NRKR=RR(1,K)IF(KR.EQ.I) C=RR(2,K)ENDDONC=C !NC=0,提取了整数部分C=C-NC !C=0.***,例如C=0.111DO J=1,NFC=C*10.0 !例如C=1.21L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位 !上的数字,这里"+0.1"是为了防止四舍五入是出现错误 C=C-LIF(L.EQ.0)THENN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DO IE=1,NEDO I=1,NDNI=ME(I,IE)DO J=1,NFLMT((I-1)*NF+J,IE)=IT(J,NI)ENDDOENDDOENDDORETURNEND!二维总刚中对角线元位置数组SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)DIMENSION MAXA(NPF),LMT(NDF,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NNMIP=I-1IG=IPDO IE=1,NEDO J=1,NDFIF(LMT(J,IE).EQ.IP) THENDO K=1,NDFIF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE)ENDDOEND IFENDDOENDDOMAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK= MAXA(NNM)-1RETURNEND!生成一维存储结构总刚度矩阵SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM)DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),& MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)CKK=0DO 10 IE=1,NETAK=0CALL 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,6DO 220 J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IF((NJ-NI).GE.0.AND.NI*NJ.GT.0) THENIJ=MAXA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIF220 CONTINUE10 CONTINUERETURNEND!生成荷载矩阵SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF)V=0PP=0DO I=1,NFDO J=1,NPDO K=1,NCFIF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J))=PF(I+1,K)ENDDOENDDODO K=1,NCFDO I=1,NPIF(I.EQ.PF(1,K))THENPP(NF*(I-1)+1)=PF(2,K)PP(NF*(I-1)+2)=PF(3,K)PP(NF*(I-1)+3)=PF(4,K)ENDIFENDDOENDDORETURNEND!对一维结构总刚度矩阵进行矩阵分解（LDLT）SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM) DIMENSION A(NWK),MAXA(NNM)IF(NN.EQ.1) RETURNDO 200 N=1,NNKN=MAXA(N)KL=KN+1KU=MAXA(N+1)-1KH=KU-KLIF(KH)304,240,210210 K=N-KHIC=0KLT=KUDO 260 J=1,KHKLT=KLT-1IC=IC+1KI=MAXA(K)ND=MAXA(K+1)-KI-1IF(ND) 260,260,270270 KK=MIN0(IC,ND)C=0.0DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L)A(KLT)=A(KLT)-C260 K=K+1240 K=NB=0.0DO 300 KK=KL,KUK=K-1KI=MAXA(K)C=A(KK)/A(KI)IF(ABS(C).LT.1.0E+07) GOTO 290WRITE(IOUT,2010) N,CSTOP290 B=B+C*A(KK)300 A(KK)=CA(KN)=A(KN)-B304 IF(A(KN)) 310,310,200310 IF(ISH.EQ.0) GOTO 320IF(A(KN).EQ.0.0) A(KN)=-1.0E-16GOTO 200320 WRITE(IOUT,2000) N,A(KN)STOP200 CONTINUERETURN2000 FORMAT(//,' Stop-stiffness matrix not positive definite',//,'no positive& pivot for equation',I4,&//,' pivot =',E20.10)2010 FORMAT(//,' Stop-strum sequence check failed + because of multiplier& growth for column &number',I4,//,' Multiplier = ',E20.8)END!回代，求得节点位移SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM)DIMENSION A(NWK),V(NN,1),MAXA(NNM)NIP=1DO IP=1,NIPDO 400 N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 400,410,410410 K=NC=0.0DO 420 KK=KL,KUK=K-1420 C=C+A(KK)*V(K,IP)V(N,IP)=V(N,IP)-C400 CONTINUEDO 480 N=1,NNK=MAXA(N)480 V(N,IP)=V(N,IP)/A(K)IF(NN.EQ.1)RETURNN=NNDO 500 L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 500,510,510510 K=NDO 520 KK=KL,KUK=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1ENDDORETURNEND!求解杆件内力、支反力和位移SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y, Z,PP,FF,SG,SM,ME,&NR,RR,NF)DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),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),PP(NPF),&SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)SG=0;SM=0;FF=0;FF2=0DO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0) THENDIST(3*(I-1)+J)=0.0ELSEIF(B.LE.NN) THENDIST(3*(I-1)+J)=FTOOL(LAB)ENDIFENDDOENDDODO IE=1,NEN1=ME(1,IE);N2=ME(2,IE)DO J=1,NFUE(J)=DIST(3*(N1-1)+J)UE(3+J)=DIST(3*(N2-1)+J)ENDDOCALL 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,NFFF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOISW=NAE(IE) AO=AE(2,ISW) SG(IE)=FE1(2)SM(IE)=FE1(2)/AODO I=1,NPFFF2(I)=FF(I)-PP(I)ENDDOENDDODO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)' 'WRITE(12,*)'****************************************'WRITE(12,*)'*********The Results of Calculation**********'WRITE(12,*)'****************************************'WRITE(12,600)WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000, I=1,NP)WRITE(12,620)WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)WRITE(12,640)WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)600 FORMAT(6X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)')610 FORMAT(1X,I4,2X,1P3E12.2)620 FORMAT(//6X,'The Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'σ(MPa)') 630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)640 FORMAT(//6X,'The Bearing Force'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')650 FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E 都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（b）所示。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
空间桁架结构程序设计(F o r t r a n)空间桁架静力分析程序及算例1、变量及数组说明2、空间桁架结构有限元分析程序源代码!主程序（读入文件，调用总计算程序，输出结果）CHARACTER IDFUT*20,OUTFUT*20WRITE(*,*) 'Input Data File name:'READ (*,*)IDFUTOPEN (11,FILE=IDFUT,STATUS='OLD')WRITE(*,*) 'Output File name:'READ (*,*)OUTFUTOPEN(12,FILE=OUTFUT,STATUS='UNKNOWN')WRITE(12,*)'*****************************************'WRITE(12,*)'* Program for Analysis of Space Trusses *'WRITE(12,*)'* School of Civil Engineering CSU *'WRITE(12,*)'* 2012.6.25 Designed By MuZhaoxiang *'WRITE(12,*)'*****************************************'WRITE(12,*)' 'WRITE(12,*)'*****************************************'WRITE(12,*)'*************The Input Data****************'WRITE(12,*)'*****************************************'WRITE(12,100)READ(11,*)NF,NP,NE,NM,NR,NCF,NDWRITE(12,110)NF,NP,NE,NM,NR,NCF,ND100 FORMAT(6X,'The General Information'/2X,'NF',5X,'NP',5X,'NE',5X,'NM',5X,'NR',& 5X,'NCF',5X,'ND')110 FORMAT(2X,I2,6I7)NPF=NF*NPNDF=ND*NFCALL ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)END!******************************************************************** !总计算程序SUBROUTINE ANALYSE(NF,NP,NE,NM,NR,NCF,ND,NPF,NDF)DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(ND,NE),IT(NF,NP),RR(ND,NR), NAE(NE),&AE(1,2),PF(4,NCF),LMT(NDF,NE),MAXA(NPF),CKK(1000),V(NPF),DIST(NPF),& PP(NPF),FF(NPF),SG(NE),SM(NE)READ(11,*)(X(I),Y(I),Z(I),I=1,NP)READ(11,*)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)READ(11,*)(RR(1,J),RR(2,J),J=1,NR)READ(11,*)(AE(1,J),J=1,2)WRITE(12,120)WRITE(12,121)(I,X(I),Y(I),Z(I),I=1,NP)WRITE(12,130)WRITE(12,131)(MM(I),ME(1,I),ME(2,I),NAE(I),I=1,NE)WRITE(12,140)WRITE(12,141)(INT(RR(1,J)),RR(2,J),J=1,NR)WRITE(12,150)WRITE(12,151)(AE(1,J),J=1,2)IF(NCF/=0)THENREAD(11,*)((PF(I,J),I=1,4),J=1,NCF)WRITE(12,160)WRITE(12,161)(INT(PF(1,J)),PF(2,J),PF(3,J),PF(4,J),J=1,NCF)ENDIF120 FORMAT(/6X,'The Information of Joints'/2x,'Joint',5X,'X',5X,'Y',5X,'Z')121 FORMAT(1X,I4,3F8.1)130 FORMAT(/6X,'The Information ofMembers'/2x,'Member',2X,'START',4X,'END',6X,'NAE')131 FORMAT(1X,I4,3I8)140 FORMAT(/6X,'The Information of SUPPORTS'/2x,'Joint',5X,'S')141 FORMAT(1X,I4,F8.3)150 FORMAT(/6X,'The Information of Sections'/4x,'E0',8X,'A0')151 FORMAT(1X,1PE8.2,F8.4)160 FORMAT(/6X,'The Loading at Joints'/2x,'Joint',5X,'FX',5X,'FY',7X,'FZ')161 FORMAT(1X,I4,3F8.2)CALL FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)CALL FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)CALL LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)CALL CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM) ISH=1CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NNM)CALL REBACK(CKK,V,MAXA,NN,NWK,NNM)CALLDISPLS(NP,NE,NPF,NM,NN,IT,V,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,NR,RR,NF) END!********************************************************************!矩阵转置子程序SUBROUTINE MAT(M,N,A,B)DIMENSION A(M,N),B(N,M)DO I=1,MDO J=1,NB(J,I)=A(I,J)END DOEND DORETURNEND!单元刚度矩阵的形成SUBROUTINE FKE(NP,NE,NM,IE,X,Y,Z,ME,NAE,AE,AKE)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/BLAKE(1,1)=CAKE(1,2)=-CAKE(2,1)=-CAKE(2,2)=CRETURNEND!单元坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6)T=0N1=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)/BLCY=(Y2-Y1)/BLCZ=(Z2-Z1)/BLT(1,1)=CX;T(2,4)=CXT(1,2)=CY;T(2,5)=CYT(1,3)=CZ;T(2,6)=CZRETURNEND!生成单元联系数组LMTSUBROUTINE FLMT(NP,NE,NN,NNM,NR,RR,ND,NF,NDF,ME,IT,LMT)DIMENSION IT(NF,NP),LMT(NDF,NE),ME(ND,NE),RR(2,NR)NN=0;NNM=0;IT=0;LMT=0N=0DO I=1,NPC=0DO K=1,NRKR=RR(1,K)IF(KR.EQ.I) C=RR(2,K)ENDDONC=C !NC=0,提取了整数部分C=C-NC !C=0.***,例如C=0.111DO J=1,NFC=C*10.0 !例如C=1.21L=C+0.1 !提取C整数部分,例如L=1,即提取了约束RR(2,K)十分位 !上的数字,这里"+0.1"是为了防止四舍五入是出现错误C=C-LIF(L.EQ.0)THENN=N+1IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNNM=NN+1DO IE=1,NEDO I=1,NDNI=ME(I,IE)DO J=1,NFLMT((I-1)*NF+J,IE)=IT(J,NI)ENDDOENDDOENDDORETURNEND!二维总刚中对角线元地址数组SUBROUTINE FMAXA(NNM,NE,LMT,MAXA,NWK,NPF,NDF)DIMENSION MAXA(NPF),LMT(NDF,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NNMIP=I-1IG=IPDO IE=1,NEDO J=1,NDFIF(LMT(J,IE).EQ.IP) THENDO K=1,NDFIF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE)ENDDOEND IFENDDOENDDOMAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK= MAXA(NNM)-1RETURNEND!生成一维存储结构总刚度矩阵SUBROUTINE CONKB(NP,NE,NM,NWK,ME,X,Y,Z,AE,NAE,LMT,MAXA,CKK,NNM) DIMENSIONCKK(NWK),X(NP),Y(NP),Z(NP),AE(2,NM),NAE(NE),LMT(6,NE),ME(2,NE),&MAXA(NNM),AK(6,2),AKE(2,2),T(2,6),TT(6,2),TAK(6,6)CKK=0DO 10 IE=1,NETAK=0CALL 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 J=1,6NI=LMT(I,IE)NJ=LMT(J,IE)IF((NJ-NI).GE.0.AND.NI*NJ.GT.0) THENIJ=MAXA(NJ)+NJ-NICKK(IJ)=CKK(IJ)+TAK(I,J)ENDIF220 CONTINUE10 CONTINUERETURNEND!生成荷载矩阵SUBROUTINE LP(V,PP,IT,PF,NN,NCF,NF,NP,NPF)DIMENSION V(NN),PP(NPF),IT(NF,NP),PF(4,NCF)V=0PP=0DO I=1,NFDO J=1,NPDO K=1,NCFIF(J.EQ.PF(1,K).AND.IT(I,J).NE.0)THENV(IT(I,J))=PF(I+1,K)ENDIFENDDOENDDOENDDODO K=1,NCFDO I=1,NPIF(I.EQ.PF(1,K))THENPP(NF*(I-1)+1)=PF(2,K)PP(NF*(I-1)+2)=PF(3,K)PP(NF*(I-1)+3)=PF(4,K)ENDIFENDDOENDDORETURNEND!对一维结构总刚度矩阵进行矩阵分解（LDLT）SUBROUTINE LDLT(A,MAXA,NN,ISH,IOUT,NWK,NNM) DIMENSION A(NWK),MAXA(NNM)IF(NN.EQ.1) RETURNDO 200 N=1,NNKN=MAXA(N)KU=MAXA(N+1)-1KH=KU-KLIF(KH)304,240,210210 K=N-KHIC=0KLT=KUDO 260 J=1,KHKLT=KLT-1IC=IC+1KI=MAXA(K)ND=MAXA(K+1)-KI-1IF(ND) 260,260,270270 KK=MIN0(IC,ND)C=0.0DO 280 L=1,KK280 C=C+A(KI+L)*A(KLT+L)A(KLT)=A(KLT)-C260 K=K+1240 K=NB=0.0DO 300 KK=KL,KUK=K-1KI=MAXA(K)C=A(KK)/A(KI)IF(ABS(C).LT.1.0E+07) GOTO 290WRITE(IOUT,2010) N,CSTOP290 B=B+C*A(KK)300 A(KK)=CA(KN)=A(KN)-B304 IF(A(KN)) 310,310,200310 IF(ISH.EQ.0) GOTO 320IF(A(KN).EQ.0.0) A(KN)=-1.0E-16GOTO 200320 WRITE(IOUT,2000) N,A(KN)STOP200 CONTINUERETURN2000 FORMAT(//,' Stop-stiffness matrix not positive definite',//,'no positive& pivot for equation',I4,&//,' pivot =',E20.10)2010 FORMAT(//,' Stop-strum sequence check failed + because of multiplier& growth for column &number',I4,//,' Multiplier = ',E20.8)END!回代，求得节点位移SUBROUTINE REBACK(A,V,MAXA,NN,NWK,NNM)DIMENSION A(NWK),V(NN,1),MAXA(NNM)NIP=1DO IP=1,NIPDO 400 N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 400,410,410410 K=NC=0.0DO 420 KK=KL,KUK=K-1420 C=C+A(KK)*V(K,IP)V(N,IP)=V(N,IP)-C400 CONTINUEDO 480 N=1,NNK=MAXA(N)480 V(N,IP)=V(N,IP)/A(K)IF(NN.EQ.1)RETURNN=NNDO 500 L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 500,510,510510 K=NDO 520 KK=KL,KUK=K-1520 V(K,IP)=V(K,IP)-A(KK)*V(N,IP)500 N=N-1ENDDORETURNEND!求解杆件内力、支反力和位移SUBROUTINE DISPLS(NP,NE,NPF,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM,ME,&NR,RR,NF)DIMENSION IT(3,NP),DIST(NPF),FTOOL(NPF),X(NP),Y(NP),Z(NP),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),PP(NPF),& SG(NE),SM(NE),FF2(NPF),RR(2,NR),FL(3*NR)SG=0;SM=0;FF=0;FF2=0DO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0) THENDIST(3*(I-1)+J)=0.0ELSEIF(B.LE.NN) THENDIST(3*(I-1)+J)=FTOOL(LAB)ENDIFENDDOENDDODO IE=1,NEN1=ME(1,IE);N2=ME(2,IE)DO J=1,NFUE(J)=DIST(3*(N1-1)+J)UE(3+J)=DIST(3*(N2-1)+J)ENDDOCALL 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,NFFF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOISW=NAE(IE)AO=AE(2,ISW)SG(IE)=FE1(2)SM(IE)=FE1(2)/AODO I=1,NPFFF2(I)=FF(I)-PP(I)ENDDOENDDODO I=1,NPDO J=1,NFLAB=IT(J,I)IF(LAB.EQ.0)THENK=K+1FL(K)=FF2(3*(I-1)+J)ENDIFENDDOENDDOWRITE(12,*)' 'WRITE(12,*)'****************************************' WRITE(12,*)'*********The Results of Calculation**********' WRITE(12,*)'****************************************'WRITE(12,600)WRITE(12,610)(I,DIST(3*I-2)*1000,DIST(3*I-1)*1000,&DIST(3*I)*1000, I=1,NP)WRITE(12,620)WRITE(12,630)(IE,SG(IE),SM(IE)/1000,IE=1,NE)WRITE(12,640)WRITE(12,650)(INT(RR(1,I)),FL(3*I-2),FL(3*I-1),FL(3*I),I=1,NR)600 FORMAT(6X,'The Joint Displacement'/2x,'Joint',6X,'X(mm)',8X,'Y(mm)',6X,'Z(mm)') 610 FORMAT(1X,I4,2X,1P3E12.2)620 FORMAT(//6X,'The Terminal Forces'/2x,'Member', 6X,'FN(kN)',6X,'σ(MPa)')630 FORMAT(3X,I4,2X,F8.2,6X,F8.2)640 FORMAT(//6X,'The Bearing Force'/2x,'Joint',8X,'X',8X,'Y',8X,'Z')650 FORMAT(2X,I4,2X,3F10.2)RETURNEND3、算例以下图所示空间桁架为例：圆形桁架穹项，其几何尺寸如图（a）所示，整体坐标系原点取在拱顶，集中荷载P作用于拱顶，各杆截面面积A和弹性模量E都相同（取E=210GPa，A=0.04m2）；各杆件及结点编号如图（b）所示。