空间桁架结构程序设计(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. 编程软件的选择：选择合适的编程软件，例如Rhino、Grasshopper等，以及相应的插件或脚本工具，为编程工作提供支持。
3. 编程语言的学习：学习一种适合新代桁架系统编程的语言，例如Python、C#等，掌握其基本语法和常用函数，为编程工作做好准备。
三、新代桁架系统编程的实践步骤1. 数据准备：收集和整理与新代桁架系统相关的数据，包括结构参数、材料属性、受力要求等，为后续的编程工作提供输入。
2. 桁架生成：通过编程软件，根据输入的数据和设计要求，生成新代桁架系统的初始结构模型，并进行初步的受力分析。
3. 结构优化：通过编程工具，对新代桁架系统进行优化设计，以提高结构的性能和效率，例如减少材料使用量、优化节点连接等。
4. 受力分析：利用编程软件进行新代桁架系统的受力分析，验证结构的稳定性和强度，确保其满足设计要求。
5. 结果展示：通过编程软件将优化后的新代桁架系统进行可视化展示，以便设计师和工程师们更好地理解和评估结构的性能。
四、新代桁架系统编程的应用案例1. 建筑结构设计：通过新代桁架系统的编程，实现建筑结构的设计和优化，提高结构的稳定性和安全性。
2. 建筑外观设计：利用新代桁架系统的编程，实现建筑外观的设计和变形效果，为建筑增添独特的艺术魅力。
结构矩阵分析与程序设计(桁架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。
空间桁架结构程序的设计(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 在局部坐标系下的刚度矩阵。
在空间桁架中，各杆方向不尽相同。
每根杆件采用各自的局部坐标系，这对于单元分析十分方便。
有限元大作业桁架课设
有限元基础大作业空间桁架结构分析学号：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本题中采用的为置大数法，即将给定的节点位移处的结构刚度矩阵中对应行的主元素换成计算机所能允许的大数。
空间桁架结构程序设计(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）所示。
桁架建筑结构设计方案
桁架建筑结构设计方案桁架结构是一种常见的建筑结构形式，它由一系列的梁和柱组成，通过形成三角形的稳定结构来承载荷载。
桁架结构具有一定的优势和特点，广泛应用于建筑设计中。
本文将介绍桁架结构设计方案，并探讨其特点和应用。
桁架结构设计方案的基本原理是利用三角形的稳定性。
通过将梁和柱组合形成不同形式的三角形结构，可以使结构更加稳定，减少材料的使用量。
桁架结构在构造上有很大的灵活性，可以根据不同的需求进行优化设计，满足不同场所的要求。
桁架结构的设计方案需要考虑以下几个方面。
首先是荷载分析。
根据建筑物的使用要求和地理条件，确定所需承载的重量和力。
结构设计师需要计算荷载的大小和方向，以确定梁柱的位置和尺寸。
其次是结构的形式和材料选择。
桁架结构可以有多种形式，包括平面桁架、空间桁架和曲面桁架等。
根据具体需求和建筑物的形状，选择相应的结构形式。
材料的选择也十分重要，需要考虑材料的强度、稳定性和耐久性等因素。
桁架结构设计方案的特点有很多。
首先是结构的轻量化。
相比于传统的混凝土结构或砖石结构，桁架结构采用金属材料或木材材料，具有更轻的重量。
这使得构造更加便捷，减少了对基础的要求，降低了建设成本。
其次是结构的坚固性。
桁架结构采用三角形的稳定结构，使得整个建筑物能够更好地抵抗外部荷载的作用，具有更好的抗震性能。
同时，桁架结构还具有可拆卸和可移动的特点，方便日后的维护和改造。
桁架结构的应用非常广泛。
在工业建筑中，桁架结构常用于机场、体育馆和仓库等大跨度建筑的设计。
由于桁架结构具有强度高、承载能力大的优势，适合于大跨度结构的设计。
此外，桁架结构还常用于桥梁、塔架和天线等工程项目的建设，能够满足大跨度结构的要求。
在特殊环境下，如地震区域或多风区域，采用桁架结构可以提高建筑物的抗震性能和风力稳定性。
总之，桁架结构设计方案是一种应用广泛的建筑结构形式。
它利用三角形的稳定性和优秀的性能，能够满足不同场所和条件下的建筑需求。
在未来的建筑设计中，桁架结构将继续发挥其独特的优势，为建筑行业做出更大的贡献。
Fortran语言编写的有限元结构程序
算例一计算简图及结果输出用平面刚架静力计算程序下图结构的内力。
各杆EA，EI相同。
已知：642EA=4.010KN,EI=1.610KN m⨯⨯∙计算简图如下：（1）输入原始数据控制参数3,5,8,7,1,2（NE,NJ,N,NW,NPJ,NPF）结点坐标集结点未知量编号0.0,0.0,0,0 0.0,4.0,1,2,3 0.0,4.0,1,2,4 4.0,4.0,5,6,7 4.0,0.0,0,0,8单元杆端结点编号及单元EA、EI 1,2,4.0E+06,1.6E+04 3,4,4.0E+06,1.6E+04 5,4,4.0E+06,1.6E+04结点荷载7.0,-15.非结点荷载1.0,2.0,2.0,-2.0,1.0,4.0,-25.0（2）输出结果NE= 3 NJ= 5 N= 8 NW= 7 NPJ= 1 NPF= 2 NODE X Y XX YY ZZ1 0.0000 0.0000 0 0 02 0.0000 4.0000 1 2 33 0.0000 4.0000 1 2 44 4.0000 4.000056 75 4.0000 0.0000 0 0 8ELEMENT NODE-I NODE-J EA EI1 12 0.400000E+07 0.160000E+052 3 4 0.400000E+07 0.160000E+053 54 0.400000E+07 0.160000E+05CODE PX-PY-PM7. -15.0000ELEMENT IND A Q1. 2. 2.0000 -18.00002. 1. 4.0000 -25.0000NODE U V CETA1 0.000000E+00 0.000000E+00 0.000000E+002 -0.221743E-02 -0.464619E-04 -0.139404E-023 -0.221743E-02 -0.464619E-04 0.357876E-024 -0.222472E-02 -0.535381E-04 -0.298554E-025 0.000000E+00 0.000000E+00 0.658499E-03ELEMENT N Q M1 N1= 46.4619 Q1= 10.7119 M1= -6.8477N2= -46.4619 Q2= 7.2881 M2= 0.00002 N1= 7.2881 Q1= 46.4619 M1= 0.0000N2= -7.2881 Q2= 53.5381 M2= 14.15233 N1= 53.5381 Q1= 7.2881 M1= 0.0000N2= -53.5381 Q2= -7.2881 M2= -29.1523算例二计算简图及结果输出用平面刚架静力计算程序下图结构的内力。
- 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）所示。