空间桁架结构程序的设计(Fortran)
维斯塔斯矩形钢管空间桁架连廊结构设计(全文)
维斯塔斯矩形钢管空间桁架连廊结构设计(全文)范本一：正文：一、引言维斯塔斯矩形钢管空间桁架连廊结构设计主要目的是为了满足建筑物中连廊结构的承重和抗震需求。
本文将详细介绍该结构的设计思路、参数计算、构件选择、荷载分析等内容。
二、设计思路该连廊结构采用矩形钢管空间桁架作为主要承重结构，以满足建筑物连廊的跨度要求。
设计思路主要包括结构形式的选择、受力分析和稳定性分析。
2.1 结构形式选择在连廊结构的设计中，考虑到跨度较大，矩形钢管空间桁架结构能够在保证结构稳定性的同时满足承重要求。
因此选择矩形钢管空间桁架结构作为主要承重结构。
2.2 受力分析在受力分析中，首先需要计算连廊结构的自重荷载。
然后考虑到连廊上可能发生的活载荷载和风荷载，进行荷载分析。
最后结合连廊的抗震设计，确定连廊结构的主要受力情况。
2.3 稳定性分析稳定性分析是为了保证连廊结构在使用过程中不发生倾覆或失稳。
需要考虑到连廊结构的刚度，通过横向稳定分析和纵向稳定分析，确定连廊结构的稳定性。
三、参数计算参数计算是在设计中必不可少的环节，包括截面尺寸、材料强度、构件节点的设计等。
3.1 截面尺寸根据连廊的跨度和荷载情况，计算连廊结构所需的最大截面尺寸。
一般情况下，矩形钢管的高度和宽度需要满足一定的宽高比要求，以保证结构的稳定性。
3.2 材料强度在材料的选择中，需要考虑到矩形钢管的强度、刚度和耐久性等。
通过计算材料的抗弯强度、抗压强度和抗拉强度，确定矩形钢管的材料强度。
3.3 构件节点设计构件节点的设计是确保连廊结构的节点连接紧固可靠、不发生脱开或错位的重要环节。
通过合理的节点设计，保证矩形钢管的连接稳定性。
四、荷载分析荷载分析是为了确定连廊结构的最大受力情况，包括自重荷载、活载荷载和风荷载。
4.1 自重荷载自重荷载主要考虑连廊结构本身的重量。
根据材料的密度和结构的截面尺寸，计算出连廊结构的自重荷载。
4.2 活载荷载活载荷载主要考虑连廊上可能承载的人员和设备等活动荷载。
桁架搭建方案
桁架搭建方案第1篇桁架搭建方案一、项目背景随着我国经济的快速发展，各类建筑工程对结构安全性、经济性和美观性的要求越来越高。
桁架结构作为一种常见的建筑结构形式，具有受力合理、节约材料、施工速度快等优点，被广泛应用于各类建筑工程中。
为确保桁架搭建过程的合法合规，提高工程质量和安全性，特制定本桁架搭建方案。
二、桁架结构选型1. 根据工程需求，结合建筑物的用途、跨度、荷载等因素，选择合适的桁架结构形式。
2. 桁架结构形式包括：平面桁架、空间桁架、空腹桁架、张弦桁架等。
3. 本项目采用平面桁架结构，桁架节点采用焊接连接，确保结构安全可靠。
三、材料选择与要求1. 钢材：应符合国家相关标准规定，具有出厂合格证及检测报告。
2. 钢材品种、规格、性能等应符合设计要求。
3. 焊接材料：应选用与母材性能匹配的焊接材料，确保焊接质量。
4. 混凝土：应符合国家现行标准规定，混凝土强度等级应符合设计要求。
四、施工工艺1. 施工准备：对施工现场进行清理，确保施工场地平整、干净。
2. 放线定位：根据设计图纸，放出桁架的轴线、边界线及预埋件位置线。
3. 钢材下料：根据设计图纸和施工工艺，进行钢材下料。
4. 钢材加工：对下料的钢材进行加工，包括切割、弯曲、焊接等。
5. 桁架拼装：在施工现场进行桁架的拼装，拼装顺序应符合施工工艺要求。
6. 焊接施工：采用二氧化碳气体保护焊进行焊接，确保焊接质量。
7. 桁架吊装：采用合适的吊装设备进行桁架的吊装，吊装过程应符合安全规范。
8. 质量检测：对桁架结构进行质量检测，包括尺寸偏差、焊缝质量等。
9. 防腐处理：对桁架结构进行防腐处理，确保结构使用寿命。
10. 结构验收：结构验收应符合国家现行标准规定，验收合格后方可进行后续施工。
五、施工安全措施1. 施工前，对施工人员进行安全技术交底，确保施工人员了解施工过程中的安全注意事项。
2. 施工现场应设置安全警示标志，提醒施工人员注意安全。
3. 施工过程中，应定期对施工设备进行检查、维护，确保设备安全运行。
平面桁架程序计算原理及程序编制
对于该平面桁架，有
杆号 L
△L
1
2
3
4
a
a
a
sqrt(2)a
v1
u1-u2
v2
(u1+ v1)/ sqrt(2)
把上述值代入应变能表达式，得到
5 sqrt(2)a (v2- u2)/ sqrt(2)
U E 2 F a v 1 2 ( u 1 u 2 ) 2 v 2 2 ( u 1 v 2 ) 2 /2 2 ( v 2 u 2 ) 2 /2 2
的内力。
5-5 结构计算简图的数据结构
完整而确切描述一个平面桁架结构的数据有三个方面：
（1）结构本体描述数据（NW, IESG, NU, X, Y, HL, HR）
（2）性质数据（EF）
（3）荷载数据（PX, PY）
➢NW为节点总数 ➢IESG杆件总数 ➢NU可动节点总数 ➢X, Y 节点坐标 ➢HL, HR 每根杆件两端节点编号 ➢EF 性质数据 ➢PX, PY 外载荷数据
FAFBN
FA和FB在x轴、y轴方向的分量分别为：
F A xF Ac o s N c o s
F A yF AsinNsin
F B xF BcosN cos
F B yF BsinN sin
把 F A x 、F B x 、 F A y 、 F B y 排列成列向量 P ，则有
FAx cos
式中 L 为杆的伸长， L 为杆的长度， N 为杆的
内力，EF / L 称为单根杆的刚度（单元刚度阵）。
注意⊿L= uB-uA，并且由平衡关系得到
FAFBNE L F(uAuB)
可把FA、FB排成列向量{N}，uA， uB，排成列向 量{u}，系数阵排成矩阵RD，即
空间结构内力位移计算
作业四编制有限元程序求解结构内力位移第一部分程序设计过程和子程序的说明本例为空间桁架结构有限元分析程序。
设计思路为：自然离散桁架结构，确定各节点自由度；为单元和节点编号，输入支撑信息、荷载信息、截面特性；运行程序求得刚度矩阵，继而求得节点位移和杆件内力。
程序所需各子程序已给出，在此只需编制主程序然后调用子程序求解即可。
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工程主体为钢筋混凝土结构，屋盖为空间管桁架结构，桁架采用稳定性比较好的倒三角形结构体系。
空间桁架设计方法
空间桁架设计方法摘要:本文将空间桁架设计方法引入到钢结构胶带机通廊的结构设计工作中，较为真实的反映了钢结构通廊的空间工作状态，为今后进行钢结构胶带机通廊的结构设计起到了一定的指导作用。
关键词:胶带机通廊;钢结构;空间桁架;结构设计Abstract: In this paper, the space truss design method is introduced to the steel structure conveyor gallery structure design work, it is a true reflection of the steel structure gallery space working condition, for future steel conveyor gallery structure design can play a guiding role.Key words: conveyor Gallery; steel structure; space truss; structure design1.引言胶带机由于具有运送散粒物料输送量大、可实现连续供料等优点，现在已经逐渐成为冶金行业运送各种散粒物料的主要设施。
随着冶金企业规模的扩大，冶金厂区建（构）筑物日益密集，加之生产工艺日趋复杂，胶带机通廊正向着大跨、超高、重载的方向发展。
而钢结构则以其轻质高强、跨度大（跨度60m左右的胶带机通廊在冶金行业的上料系统已经很常见【文献1】）、施工周期短等优点逐渐成为胶带机通廊设计的首选结构形式。
钢结构胶带机通廊主体部分由两侧承重主桁架、主桁架上弦支撑（垂直支撑和斜向支撑）、主桁架下弦支撑（垂直支撑和斜向支撑）组成，主桁架上弦支撑、主桁架下弦支撑通过节点板或螺栓将两侧承重主桁架连接成整体的空间桁架结构。
以笔者所在单位为例，钢结构胶带机通廊主体部分在结构设计阶段，空间桁架整体结构分析时将空间桁架简化为平面桁架进行内力分析：两侧承重主桁架采用PKPM系列结构分析软件中STS钢结构桁架进行建模分析；主桁架上弦支撑、主桁架下弦支撑中的垂直支撑由于分别承担屋面檩条和胶带机支腿荷载（图1），依据【文献2】按压弯构件进行复核；主桁架上弦支撑、主桁架下弦支撑中的斜向支撑依据【文献2】按拉、压杆进行复核，综合考虑构造要求进行设计。
桁架有限元程序流程（有限元课程设计）
有限单元法课程设计有限单元法是基于连续介质力学基础上发展起来的,目前使用最广泛的数值计算方法。
有限单元法解决问题的前提是各单元相邻边界的位移协调。
有限单元法解决问题的前提是各单元相邻边界的位移协调。
有限单元有限单元法将连续的求解域离散为一组有限个单元组成的组合体,由细分单元去逼近求解域,由于单元的不同连接方式和形式各异的单元形状由于单元的不同连接方式和形式各异的单元形状,,因此可以适应几何形状复杂的求解区域杂的求解区域;;第二第二,,利用每一个单元内的近似函数利用每一个单元内的近似函数((形函数形函数))来表示全求解域上待求的未知场函数待求的未知场函数,,把一个连续的无限自由度问题变成离散的有限自由度问题,只要求出单元结点的物理量只要求出单元结点的物理量,,就可以确定单元组合体上的其他未知场函数就可以确定单元组合体上的其他未知场函数,,如果选择合适的形函数选择合适的形函数,,随着网格密度的减小随着网格密度的减小,,近似解将逐步趋向精确解近似解将逐步趋向精确解;;第三第三,,有限单元法计算得到的总体刚度矩阵为稀疏带状矩阵,这样借助于电子计算机存储和计算的效率大大提高计算的效率大大提高,,便于处理大规模问题。
便于处理大规模问题。
从上述有限单元法的特性可知从上述有限单元法的特性可知,,其计算原理简单其计算原理简单,,但由于单元连接方式和单元形状的多元化元形状的多元化,,以及近似函数的选择合适与否以及近似函数的选择合适与否,,使得有限元法在针对具体问题求解时比较烦琐求解时比较烦琐,,正是基于这样的应用背景正是基于这样的应用背景,,本论文提出了一种更简单实用的单元模型—平面等效桁架单元模型。
最后最后,,编制有限元分析程序编制有限元分析程序,,将这种桁架单元模型运用于钢筋混凝土结构中型运用于钢筋混凝土结构中,,模型中混凝土采用等效桁架单元模型中混凝土采用等效桁架单元,,钢筋采用一维杆单元单元,,利用混凝土等效的应力应变关系对各种构件进行弹塑性分析,并试探性的提出了单元破坏准则。
空间桁架程序设计---精品资料
第六章 空间桁架程序设计第一节 概述一、计算模型及计算方法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 在局部坐标系下的刚度矩阵。
在空间桁架中，各杆方向不尽相同。
每根杆件采用各自的局部坐标系，这对于单元分析十分方便。
结构矩阵分析与程序设计(桁架vb代码)
'========================================'Structural Analysis Program For Plane Frame'========================================Option ExplicitPublic nn As Integer, ne As Integer, nd As Integer, ndf As IntegerPublic nf As Integer, npj As Integer, n As IntegerPublic al(50) As Double, t(4, 4) As Double, x(40) As Double, y(40) As DoublePublic jl(50) As Integer, jr(50) As Integer, e(50) As Double, a(50) As DoublePublic c(4, 4) As Double, r(120, 120) As Double, p(120) As Double, pe(120) As DoublePublic ibd(20) As Integer, ii(4) As Integer, bd(20) As Double, ff(4) As DoublePublic mj(20) As Integer, qj(20, 2) As Double, f(4) As Double, dis(4) As DoublePublic mf(30) As Integer, ind(30) As Double, aq(30) As Double, bq(30) As Double'==========='main program'===========Sub truss()Open "C:\Documents and Settings\Administrator\桌面\juzheng\桁架\tr3.2.txt" For Input As #1 Open "C:\Documents and Settings\Administrator\桌面\juzheng\fw.txt" For Output As #2Call input1Call wstiffCall loadCall boundCall gaussCall nqmClose 1Close 2End Sub'==============================='SUB-1 Read And Print Intial Data'==============================Sub input1()Dim inti As Integer, intj As Integer, i As Integer, j As Integer, k As IntegerDim dx, dy As DoublePrint #2, "Plane truss structural Analysis"Print #2, "*******************************"Print #2, "input data"Print #2, "= = = = ="Print #2,Print #2, "structural control data"Print #2, "---------------------"Print #2, "nn"; Spc(3); "nf"; Spc(3); "nd"; Spc(3); "ndf"; Spc(3); "ne"; Spc(2); "npj"; Spc(2);Input #1, nn, nf, nd, ndf, ne, npjn = 2 * (nn - nf)Print #2, nn; Spc(2); nf; Spc(2); nd; Spc(2); ndf; Spc(2); ne; Spc(2); npj; Spc(2); nPrint #2,Print #2, "Nodal coordinates"Print #2, "---------------------"Print #2, "Node "; Spc(2); "x"; Spc(5); "y"i = nnFor inti = 1 To iInput #1, inti, x(inti), y(inti)Print #2, inti; Spc(2); x(inti); Spc(3); y(inti)Next intiPrint #2,Print #2, "Element Information"Print #2, "---------------------"Print #2, "Ele.No."; Spc(4); ; "jl"; Spc(4); "jr"; Spc(6); "e"; Spc(10); "a"; Spc(6); "al"i = neFor inti = 1 To iInput #1, inti, jl(inti), jr(inti), e(inti), a(inti)Next intiFor inti = 1 To iIf jl(inti) >= jr(inti) Then StopNext intiFor inti = 1 To ij = jl(inti)k = jr(inti)dx = x(k) - x(j)dy = y(k) - y(j)al(inti) = Sqr(dx * dx + dy * dy)Print #2, Spc(3); inti; Spc(4); jl(inti); Spc(3); jr(inti); Spc(2); e(inti); Spc(2); a(inti); Spc(2); al(inti)Next intiPrint #2,k = npjIf k <> 0 ThenPrint #2, "Nodal Load"Print #2, "---------------------"Print #2, "i"; Spc(13); "mj"; Spc(3); "xd"; Spc(2); "yd"For inti = 1 To kInput #1, inti, mj(inti), qj(inti, 1), qj(inti, 2)Print #2, inti; Spc(1), mj(inti); Spc(1); qj(inti, 1); Spc(1); qj(inti, 2)Next intiEnd IfPrint #2,j = ndfIf j <> 0 ThenPrint #2, "Bonundary conditions"Print #2, "---------------------"Print #2, "i"; Spc(5); "ibd"; Spc(3); "bd"For inti = 1 To jInput #1, inti, ibd(inti), bd(inti)Print #2, inti; Spc(3); ibd(inti); Spc(3); bd(inti)Next intiEnd IfEnd Sub'======================================================== 'sub-2 Assemnble Structural Stiffness Matrix{R}'======================================================== Sub wstiff()Dim i As Integer, j As Integer, ie As Integer, k1 As Integer, k2 As IntegerFor i = 1 To nFor j = 1 To nr(i, j) = 0Next jNext iie = 1Do While ie <= neCall stiff(ie)Call locat(ie)For k1 = 1 To 4i = ii(k1)If i <= n ThenFor k2 = k1 To 4j = ii(k2)If j <= n Thenr(i, j) = r(i, j) + c(k1, k2)End IfNext k2End IfNext k1ie = ie + 1LoopFor i = 2 To nFor j = 1 To (i - 1)r(i, j) = r(j, i)Next jNext iEnd Sub'========================================================'sub-3 set up Stiffness Matrix[c]'========================================================Sub stiff(ie)Dim i As Integer, j As IntegerDim cx As Double, cy As Double, b1 As Double, b2 As Double, b3 As Double, b4 As Double Dim s1 As Double, s2 As Double, s3 As Double, s4 As Double, s5 As Double, s6 As Doublei = jl(ie)j = jr(ie)cx = (x(j) - x(i)) / al(ie)cy = (y(j) - y(i)) / al(ie)b1 = e(ie) * a(ie) / al(ie)s1 = cx * cx * b1s2 = b1 * cx * cys3 = cy * cy * b1c(1, 1) = s1c(1, 2) = s2c(1, 3) = -s1c(1, 4) = -s2c(2, 2) = s3c(2, 3) = -s2c(2, 4) = -s3c(3, 3) = s1c(3, 4) = s2c(4, 4) = s3For i = 2 To 4For j = 1 To (i - 1)c(i, j) = c(j, i)Next jNext iEnd Sub'======================================================== 'sub-4 set up element location vector{￠ò}'======================================================== Sub locat(ie)Dim i As Integer, j As Integeri = jl(ie)j = jr(ie)ii(1) = 2 * i - 1ii(2) = 2 * iii(3) = 2 * j - 1ii(4) = 2 * jEnd Sub'======================================================== 'SUB-5 Set Up Total Nodal Vector {P}'======================================================== Sub load()Dim i As Integer, j As Integer, k As Integeri = 1Do While i <= np(i) = 0#i = i + 1LoopIf npj > 0 ThenFor i = 1 To npjk = mj(i)p(2 * k - 1) = qj(i, 1)p(2 * k) = qj(i, 2)Next iEnd IfEnd Sub'======================================================== 'sub-8 set up coordinate transfer matrix[t]'======================================================== Sub trans(ie)Dim i As Integer, j As IntegerDim cx As Double, cy As Doublei = jl(ie)j = jr(ie)cx = (x(j) - x(i)) / al(ie)cy = (y(j) - y(i)) / al(ie)For i = 1 To 4For j = 1 To 4t(i, j) = 0#Next jNext iFor i = 1 To 4 Step 2t(i, i) = cxt(i, i + 1) = cyt(i + 1, i) = -cyt(i + 1, i + 1) = cxNext iEnd Sub'======================================================== 'sub-9 introduce support conditions'======================================================== Sub bound()Dim i As Integer, j As Integer, k As IntegerDim g As DoubleIf ndf <> 0 Theng = 1E+20For i = 1 To ndfk = ibd(i)r(k, k) = gp(k) = g * bd(i)Next iEnd If'======================================================== 'sub -10 solve equilibrium equations'======================================================== Sub gauss()Dim i As Integer, j As Integer, k As Integer, k1 As Integer, n1 As IntegerDim c As Doublen1 = n - 1For k = 1 To n1k1 = k + 1For i = k1 To nc = r(k, i) / r(k, k)p(i) = p(i) - p(k) * cFor j = k1 To nr(i, j) = r(i, j) - r(k, j) * cNext jNext iNext kp(n) = p(n) / r(n, n)For i = 1 To n1k = n - ik1 = k + 1For j = k1 To np(k) = p(k) - r(k, j) * p(j)Next jp(k) = p(k) / r(k, k)Next iFor k = 1 To neCall locat(k)For k1 = 1 To 4i = ii(k1)If i > n Thenp(i) = 0#End IfNext k1Next kPrint #2,Print #2, "output data"Print #2, "============"Print #2, "nadal displacement"Print #2, "-----------"Print #2, "node no.", Spc(5); "u", Spc(18); "v"For i = 1 To nnPrint #2, i, p(2 * i - 1), p(2 * i)NextEnd Sub'======================================================== 'sub-11 Calculate Member-end Force Of Elements'========================================================Sub nqm()Dim i As Integer, j As Integer, ie As Integer, k As IntegerDim fe(4) As Double, f1(4) As Double, fg As DoublePrint #2, "force and stress of Elements :"Print #2, "= = = = = = = = = = = = = = ="Print #2,Print #2, "Eie No.", Spc(6); "Force-N", Spc(20); "Stress-N"Print #2, "-----------------------------------"ie = 1Do While ie <= neCall trans(ie)Call stiff(ie)Call locat(ie)For i = 1 To 4dis(i) = 0#j = ii(i)If j <= n Thendis(i) = p(j)End IfNext iFor i = 1 To 4fe(i) = 0#For j = 1 To 4fe(i) = fe(i) + c(i, j) * dis(j)Next jNext iFor i = 1 To 4f(i) = 0#For j = 1 To 4f(i) = f(i) + t(i, j) * fe(j)Next jNext ifg = f(3) / a(ie)Print #2, ie, Spc(5); f(3), Spc(5); fgie = ie + 1LoopEnd Sub。
桁架结构设计步骤
桁架结构设计步骤桁架结构设计步骤如下：第一步：确定基本设计参数设计的基本参数包括板的跨度和厚度、两阶段的板支撑、钢筋类型、混凝土强度等级和使用荷载。
第二步：钢桁架楼承板长度的确定根据工程实际情况，楼承板的长度可以是一跨，也可以是多跨之和（1）钢桁架楼板的长度应为200mm的倍数，特殊情况下，长度可为100mm的倍数。
（2）楼承板的长度应为多跨之和的连续板。
（3）楼承板的长度不宜大于20m，理论上钢桁架楼承板可以加工成无限长，但实际上考虑到楼承板的运输系数，最大长度不应超过17.5米，否则很难找到运输工具。
部分项目与承重板之间不允许有严格的拼接要求。
此时，需要现场处理。
第三步：根据使用阶段计算，初步选定钢桁架楼承板的类型钢桁架楼承板设计包括四个部分：桁架构件设计、底模设计、桁架构件连接节点设计、桁架与底模连接节点设计。
其中，连接节点的强度由结构保证，无需验算。
底模设计成型，满足应力要求。
因此，设计者只需设计桁架构件就可以选择钢桁架楼承板的类型。
第四步：当没有临时支撑时，应检查表或检查施工阶段，调整地板承重板的类型，以满足应力要求。
第五步：确定支座附加钢筋的数量当钢桁架连续时，使用阶段计算的支座负筋面积减去钢桁架上弦杆截面面积，即为支座的附加配筋量；当钢桁架在支座处不连续时，支座负筋在使用阶段计算的截面面积为支座的附加钢筋用量。
不同类型的钢筋应更换为等强度带。
第六步：楼层结构图楼层结构图包括平面布置图和节点详图。
平面布置图包括：钢筋桁架楼承板、支座负筋、孔边及柱边附加钢筋、分布钢筋、柱边及混凝土墙边支撑等，同时，施工中临时支撑的布置必须在图纸中明确。
第七步：其他注意事项楼板可设计为单向板或双向板。
钢桁架楼承板在施工阶段均为单向板。
无临时支撑时，施工阶段所需钢筋一般大于使用阶段按单向板计算的钢筋，故楼板应按单向板设计。
当因具体工程条件需要设计双向板时，为节约钢材，施工阶段应沿垂直于桁架方向设置临时支撑。
空间桁架结构程序的设计(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）所示。
有限元大作业桁架课设
有限元基础大作业空间桁架结构分析学号：010810726学生姓名：同组成员：完成时间：2013年11月20日南京航空航天大学航空宇航学院二〇一一年十月八日一、问题概述该模型为一空间桁架结构，由36根截面积为10mm²，长度分别为30、40、50cm的杆组成，共有12个节点，水平面均为由30cm长杆组成的等边三角形，底部三个节点均为固定铰链约束。
杆材料相同，弹性模量为7 E+6N/cm²，µ=0.3。
在结构的5个节点上分别作用了大小均为1000N的外载荷，求相应的各干轴力。
二、模型简化（一）建立数学模型建立如下图所示的空间坐标系，以及相应的节点编号。
（二）建立有限元模型将结构离散化，根据空间桁架的受力特点将它离散为“受轴力的空间杆单元”。
所取的整体坐标、节点和单元的编号，计算所需的原始数据列表如下：单元号节点坐标lcmAcm²EN/ cm²编号i-j Xicm YicmZicmXjcmYjcmZjcm○11-4 -15 0 0 -15 0 40 40 0.1 7E+6 ○21-5 -15 0 0 0 25.98 40 50 0.1 7E+6 ○31-6 -15 0 0 15 0 40 50 0.1 7E+6 ○42-4 0 25.98 0 -15 0 40 50 0.1 7E+6 ○52-5 0 25.98 0 0 25.98 40 40 0.1 7E+6 ○62-6 0 25.98 0 15 0 40 50 0.1 7E+6 ○73-4 15 0 0 -15 0 40 50 0.1 7E+6 ○83-5 15 0 0 0 25.98 40 50 0.1 7E+6 ○93-6 15 0 0 15 0 40 40 0.1 7E+6 ○104-5 -15 0 40 0 25.98 40 30 0.1 7E+6 ○115-6 0 25.98 40 15 0 40 30 0.1 7E+6 ○126-4 15 0 40 -15 0 40 30 0.1 7E+6 ○134-7 -15 0 40 -15 0 80 40 0.1 7E+6 ○144-8 -15 0 40 0 25.98 80 50 0.1 7E+6 ○154-9 -15 0 40 15 0 80 50 0.1 7E+6 ○165-7 0 25.98 0 -15 0 80 50 0.1 7E+6 ○175-8 0 25.98 0 0 25.98 80 40 0.1 7E+6 ○185-9 0 25.98 40 15 0 80 50 0.1 7E+6 ○196-7 15 0 40 -15 0 80 50 0.1 7E+6 ○206-8 15 0 40 0 25.98 80 50 0.1 7E+6 ○216-9 15 0 40 15 0 80 40 0.1 7E+6 ○227-8 -15 0 80 0 25.98 80 30 0.1 7E+6 ○238-9 0 25.98 80 15 0 80 30 0.1 7E+6 ○249-7 15 0 80 -15 0 80 30 0.1 7E+6 ○257-10 -15 0 80 -15 0 120 40 0.1 7E+6 ○267-11 -15 0 80 0 25.98 120 50 0.1 7E+6○277-12 -15 0 80 15 0 120 50 0.1 7E+6 ○288-10 0 25.98 80 -15 0 120 50 0.1 7E+6 ○298-11 0 25.98 80 0 25.98 120 40 0.1 7E+6 ○308-12 0 25.98 80 15 0 120 50 0.1 7E+6 ○319-10 15 0 80 -15 0 120 50 0.1 7E+6 ○329-11 15 0 80 0 25.98 120 50 0.1 7E+6 ○339-12 15 0 80 15 0 120 40 0.1 7E+6 ○3410-11 -15 0 120 0 25.98 120 30 0.1 7E+6 ○3511-12 0 25.98 120 15 0 120 30 0.1 7E+6 ○3612-10 15 0 120 -15 0 120 30 0.1 7E+6（三）受载与约束情况本题所受的节点外载荷为：X7=300.00 Y7=429.60 Z7=-800.00X9=-300.00 Y9=429.60 Z9=-800.00Z10=-1000 Z11=-1000 Z12=-1000其余节点外力为0本题中采用的为置大数法，即将给定的节点位移处的结构刚度矩阵中对应行的主元素换成计算机所能允许的大数。
空间桁架
#include<stdio.h>#include<math.h>main(){inti,j,w,k,kk,v,s,r,n,m,unit[5][4]={1,1,4,1,2,2,4,1,3,3,4,1},bind[5][4]={1,0,0,0,2,0,0,0,3,0,0,0},ii[3]; float xi,yi,zi,xj,yj,zj,L,E,A,EAL,op,LEX,MEX,NEX;float point[4][3]={-4,0,0,7,-3,0,5,3,0,0,0,6},material[4][2]={2.1e9,0.0025},P[5][4]={4,0,0,-50}; float disXi,disXj,disY i,disYj,disZi,disZj,N;float unitarray[2][2][9],allarray[5*3][5*3];for(i=0;i<3*5;i++)for(w=0;w<3*5;w++)allarray[i][w]=0.0;for(k=0;k<5;k++){i=unit[k][1];j=unit[k][2];ii[0]=i;ii[1]=j;xi=point[i-1][0];yi=point[i-1][1];zi=point[i-1][2];xj=point[j-1][0];yj=point[j-1][1];zj=point[j-1][2];L=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi)+(zj-zi)*(zj-zi));E=material[unit[k][3]-1][0];A=material[unit[k][3]-1][1];EAL=E*A/L ;LEX=(xj-xi)/L;MEX=(yj-yi)/L;NEX=(zj-zi)/L;for(r=0;r<2;r++)for(s=0;s<2;s++){if(r==s){unitarray[r][s][0]=EAL*LEX*LEX;unitarray[r][s][1]=EAL*LEX*MEX;unitarray[r][s][2]=EAL*LEX*NEX;unitarray[r][s][3]=EAL*LEX*MEX;unitarray[r][s][4]=EAL*MEX*MEX;unitarray[r][s][5]=EAL*MEX*NEX;unitarray[r][s][6]=EAL*LEX*NEX;unitarray[r][s][7]=EAL*MEX*NEX;unitarray[r][s][8]=EAL*NEX*NEX;}{unitarray[r][s][0]=-EAL*LEX*LEX;unitarray[r][s][1]=-EAL*LEX*MEX;unitarray[r][s][2]=-EAL*LEX*NEX;unitarray[r][s][3]=-EAL*LEX*MEX;unitarray[r][s][4]=-EAL*MEX*MEX;unitarray[r][s][5]=-EAL*MEX*NEX;unitarray[r][s][6]=-EAL*LEX*NEX;unitarray[r][s][7]=-EAL*MEX*NEX;unitarray[r][s][8]=-EAL*NEX*NEX;}}for(v=0;v<2;v++)for(w=0;w<2;w++){allarray[3*ii[v]-3][3*ii[w]-3]+=unitarray[v][w][0]; allarray[3*ii[v]-3][3*ii[w]-2]+=unitarray[v][w][1]; allarray[3*ii[v]-3][3*ii[w]-1]+=unitarray[v][w][2]; allarray[3*ii[v]-2][3*ii[w]-3]+=unitarray[v][w][3]; allarray[3*ii[v]-2][3*ii[w]-2]+=unitarray[v][w][4]; allarray[3*ii[v]-2][3*ii[w]-1]+=unitarray[v][w][5]; allarray[3*ii[v]-1][3*ii[w]-3]+=unitarray[v][w][6]; allarray[3*ii[v]-1][3*ii[w]-2]+=unitarray[v][w][7]; allarray[3*ii[v]-1][3*ii[w]-1]+=unitarray[v][w][8];}}printf("总刚矩阵为：\n");for(r=0;r<3*4;r++){ for(s=0;s<3*4;s++)printf("%8.1g ",allarray[r][s]);printf("\n"); }printf("\n");for(i=0;i<3;i++){ if(bind[i][1]==0){for(j=0;j<3*4;j++){allarray[3*bind[i][0]-3][j]=0.0;allarray[j][3*bind[i][0]-3]=0.0;allarray[3*bind[i][0]-3][3*bind[i][0]-3]=1.0;}if(bind[i][2]==0){ for(j=0;j<3*4;j++){allarray[3*bind[i][0]-2][j]=0.0;allarray[j][3*bind[i][0]-2]=0.0;}allarray[3*bind[i][0]-2][3*bind[i][0]-2]=1.0;}if(bind[i][3]==0){for(j=0;j<3*4;j++){allarray[3*bind[i][0]-1][j]=0.0;allarray[j][3*bind[i][0]-1]=0.0;}allarray[3*bind[i][0]-1][3*bind[i][0]-1]=1.0;}}for(i=0;i<3*4;i++)allarray[i][3*4]=0.0;{allarray[3*(int)P[0][0]-3][3*4]=P[0][1];allarray[3*(int)P[0][0]-2][3*4]=P[0][2];allarray[3*(int)P[0][0]-1][3*4]=P[0][3];}n=3*4;for(k=0;k<n;k++){ op=allarray[k][k];for(m=k;m<=n;m++)allarray[k][m]=allarray[k][m]/op;for(i=k+1;i<n;i++){op=allarray[i][k];for(m=k;m<=n;m++)allarray[i][m]=allarray[i][m]-allarray[k][m]*op;}}for(k=n-1;k>=0;k--){for(i=k-1;i>=0;i--){op=allarray[i][k];for(m=n;m>=0;m--)allarray[i][m]=allarray[i][m]-allarray[k][m]*op;}}for(i=0;i<n;i++){ for(j=0;j<=n;j++)printf("%8.2g ",allarray[i][j]);printf("\n");}for(kk=0;kk<3;kk++){i=unit[kk][1];j=unit[kk][2];xi=point[i-1][0];yi=point[i-1][1];zi=point[i-1][2];xj=point[j-1][0];yj=point[j-1][1];zj=point[j-1][2];L=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi)+(zj-zi)*(zj-zi));E=material[unit[kk][3]-1][0];A=material[unit[kk][3]-1][1];EAL=E*A/L ;LEX=(xj-xi)/L;MEX=(yj-yi)/L;NEX=(zj-zi)/L;disXi=allarray[3*i-3][3*4];disYi=allarray[3*i-2][3*4];disZi=allarray[3*i-1][3*4];disXj=allarray[3*j-3][3*4];disYj=allarray[3*j-2][3*4];disZj=allarray[3*j-1][3*4];N=-EAL*(LEX*(disXi-disXj)+MEX*(disYi-disYj)+NEX*(disZi-disZj));printf("the unit=%d i=%d j=%d\n",unit[kk][0],i,j);printf("disXi=%f disY i=%f disZi=%f disXj=%f disYj=%f disZj=%f\n",disXi,disYi,disZi,disXj,disYj,disZj); printf("the inner force=%f\n",N);}}。
空间桁架结构程序设计(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）所示。
空间桁架结构程序设计(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 分析和优化设计:接下来是对初始几何模型进行分析和优化。
- 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,'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）所示。