空间桁架结构程序设计(Fortran)
维斯塔斯矩形钢管空间桁架连廊结构设计(全文)
维斯塔斯矩形钢管空间桁架连廊结构设计(全文)范本一：正文：一、引言维斯塔斯矩形钢管空间桁架连廊结构设计主要目的是为了满足建筑物中连廊结构的承重和抗震需求。
本文将详细介绍该结构的设计思路、参数计算、构件选择、荷载分析等内容。
二、设计思路该连廊结构采用矩形钢管空间桁架作为主要承重结构，以满足建筑物连廊的跨度要求。
设计思路主要包括结构形式的选择、受力分析和稳定性分析。
2.1 结构形式选择在连廊结构的设计中，考虑到跨度较大，矩形钢管空间桁架结构能够在保证结构稳定性的同时满足承重要求。
因此选择矩形钢管空间桁架结构作为主要承重结构。
2.2 受力分析在受力分析中，首先需要计算连廊结构的自重荷载。
然后考虑到连廊上可能发生的活载荷载和风荷载，进行荷载分析。
最后结合连廊的抗震设计，确定连廊结构的主要受力情况。
2.3 稳定性分析稳定性分析是为了保证连廊结构在使用过程中不发生倾覆或失稳。
需要考虑到连廊结构的刚度，通过横向稳定分析和纵向稳定分析，确定连廊结构的稳定性。
三、参数计算参数计算是在设计中必不可少的环节，包括截面尺寸、材料强度、构件节点的设计等。
3.1 截面尺寸根据连廊的跨度和荷载情况，计算连廊结构所需的最大截面尺寸。
一般情况下，矩形钢管的高度和宽度需要满足一定的宽高比要求，以保证结构的稳定性。
3.2 材料强度在材料的选择中，需要考虑到矩形钢管的强度、刚度和耐久性等。
通过计算材料的抗弯强度、抗压强度和抗拉强度，确定矩形钢管的材料强度。
3.3 构件节点设计构件节点的设计是确保连廊结构的节点连接紧固可靠、不发生脱开或错位的重要环节。
通过合理的节点设计，保证矩形钢管的连接稳定性。
四、荷载分析荷载分析是为了确定连廊结构的最大受力情况，包括自重荷载、活载荷载和风荷载。
4.1 自重荷载自重荷载主要考虑连廊结构本身的重量。
根据材料的密度和结构的截面尺寸，计算出连廊结构的自重荷载。
4.2 活载荷载活载荷载主要考虑连廊上可能承载的人员和设备等活动荷载。
MatLab与Fortran混合编程实现结构优化和可靠性分析
（ ａｕｙｏ ｏｓｕｔｎＥ ｇｎｅｎ Ｄ ｌ ｎＵ ｉｒｔ ｏ Ｔｃｎｌ ｙ ａ ａ ，Ｌ ｏｉ １０４ ｈｎ ） Ｆｃｌ ｔ ｆＣｎｔ ｃｏ ｎｉ ｒ ｇ， ａｉ ｎ ｅ ｉ ｅｈｏ ｇ ，Ｄ ｌ ｎ ｉ ｎ ｇ １６２ ，Ｃ ｉ ｒ ｉ ｅｉ ａ ｖ ｓｙｆ ｏ ｉ ａ ｎ ａ
２００ 第 第 ３ １卷 ６月 １２年 翅
Ｊｕ ｌ ｆ ｔ ｅｏｒｅ ａｄＡ ｃｉｃ 珀 ｏｍａ ０ ｗａ水利与建筑工程学报 ｅ Ｒ ｓＩｃｓ ｎ ｒ ｔ 】 ｒ ｌ ｈｅｕ
Ｖ １１ ｏ３ ｏ．０Ｎ ．
Ｊ ｎ ．， ２ ｕ ２０ １
Ｍａ ａ ｔ ｂ与 Ｆ ｒ ａ Ｌ ｏｔ ｎ混 合 编 程 实 现 ｒ 结 构 优 化 和 可 靠 性 分 析
ｓｒ ｔ ｕｃｔ ｅ ｕｒ ｓ．
Ｋ ｅ ｗｏ ｄ ｙ ｒ ｓ：Ｍ ａ Ｌａ ａ ｕ ｇ ｔ ｂ ｌｎｇ ａ ｅ；Ｆｏ ｔａ ａ ｇ ａ ｅ；ｍ ｉｉ ｏ ｒ ｍ ｍ ｉ ｇ； ｓｒ ｃ ｕ ａ ｐ ｉ ｉａ ｉ ｎ；ｓｒ ｔ ａ ｅｉ ｒｒ ｎ ｌ ｎ ｕｇ ｘ ｎｇ ｐｒ ｇ ａ ｎ ｔ ｕ ｔ ｒ ｌｏ ｔｍ ｚ ｔｏ ｔｕｃｕｒ ｌｒ ｌ－
ａ ｉｉ ｎａｙ ｉ ｂ ｌ ｙ ａ ｌ ｓｓ ｔ
Ｍ ｔａ 化 工 具 箱 针 对 各 种 优 化 问题 给 出 了 ａＬｂ优 完整 的解 决方 案 ｌ＿ ， 内容 涵盖 了包括最 值 问题 、 １ 其 ｌ
Ｆｒａ 有 限元 程序成 功 连 接 ； ｏｒ ｔｎ 然后 在 Ｍａ ａ ｔ ｂ环 境 下 ｌ
赵 秀 丽 ， 平 易
（ 大连理工大学 建设工程学部 ， 辽宁 大连 １６２ ） １０４
摘 要 ： 为有效利用 Ｍａ ａ ｔ ｂ和现有 Ｆ ￣ａ 有 限元程 序 ， Ｌ ｏｒ ｎ 采用 Ｍａ ａ 与 Ｆｒａ 合编程 。通过 Ｍａ ａ 程 ｔｂ Ｌ ｏｒ ｔ ｔｎ 昆 ｔｂ Ｌ 序接 口， Ｍａ ａ 与 Ｆ￣ａ 有 限元 程序成功连接实现平 面桁架结 构 的静 力分析 ； 将 ｔｂ Ｌ ｏｒ ｎ 然后 在 Ｍ ｔａ ａ￣ ｌｂ环境下 运用优 化工 具箱 和统计工具箱编程实现 了平面桁架 的结 构优 化设 计和可靠指标的计算。 关键词 ： ａ ａ 语言 ；ｏ ｒ 语 言 ； Ｍｔ ｂ Ｌ Ｆ￣ａ ｎ 混合编程 ； 结构优化 ； 构可靠性 分析 结
桁架有限元程序流程(有限元课程设计)
有限单元法课程设计有限单元法是基于连续介质力学基础上发展起来的,目前使用最广泛的数值计算方法。
有限单元法解决问题的前提是各单元相邻边界的位移协调。
有限单元法将连续的求解域离散为一组有限个单元组成的组合体,由细分单元去逼近求解域,由于单元的不同连接方式和形式各异的单元形状,因此可以适应几何形状复杂的求解区域;第二,利用每一个单元内的近似函数(形函数)来表示全求解域上待求的未知场函数,把一个连续的无限自由度问题变成离散的有限自由度问题,只要求出单元结点的物理量,就可以确定单元组合体上的其他未知场函数,如果选择合适的形函数,随着网格密度的减小,近似解将逐步趋向精确解;第三,有限单元法计算得到的总体刚度矩阵为稀疏带状矩阵,这样借助于电子计算机存储和计算的效率大大提高,便于处理大规模问题。
从上述有限单元法的特性可知,其计算原理简单,但由于单元连接方式和单元形状的多元化,以及近似函数的选择合适与否,使得有限元法在针对具体问题求解时比较烦琐,正是基于这样的应用背景,本论文提出了一种更简单实用的单元模型—平面等效桁架单元模型。
最后,编制有限元分析程序,将这种桁架单元模型运用于钢筋混凝土结构中,模型中混凝土采用等效桁架单元,钢筋采用一维杆单元,利用混凝土等效的应力应变关系对各种构件进行弹塑性分析,并试探性的提出了单元破坏准则。
用本文方法和商用有限元分析软件Ansys9.0的计算结果进行比较,经验证用本文模型在保证同等工程精度的条件下,是一种简单可行的方法。
关键词：有限单元法；平面桁架；形函数；刚度矩阵；有限元分析软件一、桁架有限元程序流程 (1)1、子程序说明： (1)2、平面桁架内力计算的标识符 (2)BH 二维数组，用于存放单元截面尺寸 (2)NRES 二维数组，用于存放约束的位移值 (2)JP 二维数组，用于存放节点的荷载值 (2)ESTIF 四维数组，用于存放整体坐标系下的单元刚度矩阵 (2)二、数据准备 (3)三、平面桁架内力计算程序 (4)参考文献 (11)设计题目：平面桁架程序计算单元应力设计题目：如图所示桁架，已知杆件材料的弹性模量E=2.1 104KN/cm2,杆件截面高度H=12cm，截面的宽度为b=3cm，不计各杆的自重，求在荷载作用下，各杆的轴力。
弹性力学与有限元分析第二章-平面桁架有限元分析及程序设计
x
由单元①的刚度方程：
Fj
①
k
① ji
i
①
k
① jj
j
①
k
① ji
2
k
① jj
1
由单元③的刚度方程：
Fj
③
k
③ ji
i
③
k
③ jj
j
③
k
③ ji
3
k
③ jj
1
§2.3 结点平衡与整体刚度矩阵的集成
代入结点1的平衡条件:
k
l
xi
)
(dx j
dxi
)
(
yj
l
yi )
(dy j
dyi )
(dx j dxi ) (dy j dyi )
cos sin
由于杆件的变形产生位移：
ui dxi vi dyi
u j dxj v j dy j
因此，杆件应变为：
dl l
l
(ui
uj)
l
(vi
vj)
杆件轴力为：
(2k1 k2 )v4 P
结构的整体刚度系数
v4
P 2k1
k2
12 3
l2 l1 l1
4 P
N1
N1y
cos
k1v4
cos
k1P
(2k1 k2 ) cos
N2
k2v4
k2P 2k1 k2
位移法求解超静定结构。
§2.1 平面桁架单元的离散
结构的离散化：尽量将结构离散成数量最少的等截面直 杆单元
kki③ ③jii
ki③j
k
③ jj
3 3 3 3
§2.3 结点平衡与整体刚度矩阵的集成
空间结构内力位移计算
作业四编制有限元程序求解结构内力位移第一部分程序设计过程和子程序的说明本例为空间桁架结构有限元分析程序。
设计思路为：自然离散桁架结构，确定各节点自由度；为单元和节点编号，输入支撑信息、荷载信息、截面特性；运行程序求得刚度矩阵，继而求得节点位移和杆件内力。
程序所需各子程序已给出，在此只需编制主程序然后调用子程序求解即可。
1、主要变量及数组说明程序中要设置许多变量和数组来存放各种数据。
在本程序中，变量及数组名称选用习惯中常用的表示方法，同时遵从FORTRAN90语言的隐含规则，即由字母I-N开头的均为整型，否则为实型。
程序中的变量及数组说明详见附录源程序变量和数组说明。
2、内力计算及检算程序形成一维存储总刚子程序CONKB、解线性方程组子程序LDLTREBACK、求节点位移及单元内力子程序DISPLS。
下面是子程序的介绍。
☆SUBROUTINE FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT)功能：形成IT以及LMT子程序：传入参数：NP,NE ,NR,ME,RR传出参数：IT,LMT,NN,NN1☆FMAXA(NN1,NE,LMT,MAXA,NWK,NP)功能：形成MAXA数组传入参数：NN1,NE,LMT,NP传出参数：MAXA,NWK☆SUBROUTINECONKB(NP,NE,NWK,ME,X,Y,Z,AE,LMT,MAXA,CKK,NN1)功能：①根据输入的杆件编号、节点位置、杆件位置信息及截面信息，形成杆件在局部坐标系下刚度矩阵的子程序：5SUBROUTINE FKE(IE,NP,NE,X,Y,Z,ME,AE,AKE)传入参数：NF,NP,NE,NR, ME,X,Y,Z,AE传出参数：AKE②由局部坐标系向总体坐标系转换的子程序：SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)传入参数：IE,NP,NE,X,Y,Z,ME传出参数：T③矩阵转置子程序：SUBROUTINE MAT(M,N,A,B)传入参数：M,N,A传出参数：B④矩阵乘法子程序：SUBROUTINE MUL(A,B,M,N,L,AB)传入参数：A,B,M,N,L传出参数：AB传入参数：NF,NP,NE,NM,NR,ME,X,Y,Z,AE,NAE,NN1,JS传出参数：IT,MAXA,CKK，NWK☆SUBROUTINELDLTREBACK(PF,CKK,V,MAXA,NN,NWK,NN1,ISH,IOUT,NCF,NP) 功能：①将一维存储的结构刚度矩阵进行LDLT分解子程序SUBROUTINE LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1)传入参数：CKK,MAXA,NN,ISH,IOUT,NWK,NN1传出参数：MAXA,CKK②回代求解得节点位移子程序6SUBROUTINEREBACK(PF,CKK,V,MAXA,NN,NWK,NN1,NCF,NR)传入参数：PF,CKK,MAXA,NN,NWK,NN1,NCF,NR传出参数：V☆SUBROUTINEDISPLS(NP,NE,NM,NN,IT,FTOOL,DIST,AE,NAE,X,Y,Z,PP,FF,SG,SM) 功能：求单元内力、节点位移及约束反力传入参数：NP,NE,NM,NN,IT,FTOOL,AE,NAE,X,Y,Z传出参数：DIST,SG,SM,FF第二部分程序使用说明1、在使用所述桁架内力计算程序具体步骤如下：（1）绘制结构计算简图，对杆件和节点进行编号；（2）建立整体坐标系，确定各个节点的坐标；（3）填写输入数据表；（4）建立输入数据文件，文件扩展名为“.txt”；（5）运行程序，按屏幕提示语句进行操作；（6）程序运行完毕，查看输出数据文件结果。
空间管桁架结构设计
空间管桁架结构设计王芬;周云平【摘要】空间管桁架作为一种新型的空间结构体系，具有刚度大、用钢量小、施工方便、经济环保性能优等特点，广泛应用于大型输变电站、体育馆、车站等大跨空间结构中。
通过对管桁架结构布置、稳定性和抗震性能的分析，讨论结构外观体型布置对大跨度结构内力及变形的影响，为同类结构设计提供参考和借鉴。
%As a new spatial structure system,the structure of the spatial tube truss has the features of great rigidity,small in amount of steel,convenient for constructing,economical and environmental in performance and therefore it is widely used in substations,stadiums and railway stations. Based on an analysis of the structural arrangement,stability and seismic performance the tube truss,this paper discusses the effects of the outer structural arrangement on the internal forces and deformations of the large-span structure in a bid to provide useful reference for the design of similar structures.【期刊名称】《电网与清洁能源》【年(卷),期】2015(000)007【总页数】5页(P123-127)【关键词】空间管桁架;结构布置;稳定性能;抗震性能;SAP2000设计【作者】王芬;周云平【作者单位】中煤西安设计工程有限责任公司，陕西西安 710054;西安理工大学土木建筑工程学院，陕西西安 710048【正文语种】中文【中图分类】TU323.4工程主体为钢筋混凝土结构，屋盖为空间管桁架结构，桁架采用稳定性比较好的倒三角形结构体系。
空间框架静力计算源程序(程序结构力学大作业)
主要为一些工具函数 变量定义，单元属性分析 矩阵求解 整体控制模块
1：Lxz_Tools.f90 module Lxz_Tools implicit none integer (kind(1)),parameter ::ikind=(kind(1)) integer (kind(1)),parameter ::rkind=(kind(0.D0)) real (rkind), parameter :: Zero=0.D0,One=1.D0,Two=2.D0,Three=3.D0, & & Four=4.D0,Five=5.D0,Six=6.D0,Seven=7.D0,Eight=8.D0,Nine=9.D0, & & Ten=10.D0 contains function matinv(A) result (B) real(rkind) ,intent (in)::A(:,:) !real(rkind) , allocatable::B(:,:) real(rkind) , pointer::B(:,:) integer(ikind):: N,I,J,K real(rkind)::D,T real(rkind), allocatable::IS(:),JS(:) N=size(A,dim=2) allocate(B(N,N)) allocate(IS(N));allocate(JS(N)) B=A do K=1,N D=0.0D0 do I=K,N do J=K,N if(abs(B(I,J))>D) then D=abs(B(I,J)) IS(K)=I JS(K)=J
FORTRAN语言程序设计
n 是一个十进制数字(0~9)序列。
Kind值为：1、2、4、8之一。分别表示1、2、4、8个字节个数。
FORTRAN 90/95标准中整型常量的范围没有明确规定。
例如：122、0、-36、559_2
均为合法的整型常量
100.0、100，000、100 0、1002 均为非法的整型常量
5.6.2.2 实型常量
5.5 程序的书写格式 （1）固定格式：将一个语句行分为若干个区域，如下图所示
：
5.6 数据类型
5.6.1 基本概念 1.内部数据类型 FORTRAN语言将内部数据划分为以下类型： 整型 实型 算术型 数据类型 复型 逻辑型 字符型 2.种别 一个数据在内存中均占有一定字节个数的存储单元。上述每类数据都有 其不同的种别（即Kind）特性，即上述每类数据根据其种别特性（即 Kind值）的不同分别拥有不同字节个数的存储单元。 3.数据对象 1）常量：程序运行时，其值不能改变的量，称为常量。如:5,3等等。 2）变量：程序运行时，其值可以改变的量，如:变量a, a是一个存储单元
第四章 FORTRAN语言开发环境
详见教材：3.3 FORTRAN语言 开发环境
第五章 FORTRAN 语言基础知识
5.1 字符集
FORTRAN语言允许使用的字符集为：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
空间桁架静力分析程序及算例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）所示。