空间桁架结构程序设计(Fortran)
- 格式:doc
- 大小:328.17 KB
- 文档页数:16
维斯塔斯矩形钢管空间桁架连廊结构设计(全文)范本一:正文:一、引言维斯塔斯矩形钢管空间桁架连廊结构设计主要目的是为了满足建筑物中连廊结构的承重和抗震需求。
本文将详细介绍该结构的设计思路、参数计算、构件选择、荷载分析等内容。
二、设计思路该连廊结构采用矩形钢管空间桁架作为主要承重结构,以满足建筑物连廊的跨度要求。
设计思路主要包括结构形式的选择、受力分析和稳定性分析。
2.1 结构形式选择在连廊结构的设计中,考虑到跨度较大,矩形钢管空间桁架结构能够在保证结构稳定性的同时满足承重要求。
因此选择矩形钢管空间桁架结构作为主要承重结构。
2.2 受力分析在受力分析中,首先需要计算连廊结构的自重荷载。
然后考虑到连廊上可能发生的活载荷载和风荷载,进行荷载分析。
最后结合连廊的抗震设计,确定连廊结构的主要受力情况。
2.3 稳定性分析稳定性分析是为了保证连廊结构在使用过程中不发生倾覆或失稳。
需要考虑到连廊结构的刚度,通过横向稳定分析和纵向稳定分析,确定连廊结构的稳定性。
三、参数计算参数计算是在设计中必不可少的环节,包括截面尺寸、材料强度、构件节点的设计等。
3.1 截面尺寸根据连廊的跨度和荷载情况,计算连廊结构所需的最大截面尺寸。
一般情况下,矩形钢管的高度和宽度需要满足一定的宽高比要求,以保证结构的稳定性。
3.2 材料强度在材料的选择中,需要考虑到矩形钢管的强度、刚度和耐久性等。
通过计算材料的抗弯强度、抗压强度和抗拉强度,确定矩形钢管的材料强度。
3.3 构件节点设计构件节点的设计是确保连廊结构的节点连接紧固可靠、不发生脱开或错位的重要环节。
通过合理的节点设计,保证矩形钢管的连接稳定性。
四、荷载分析荷载分析是为了确定连廊结构的最大受力情况,包括自重荷载、活载荷载和风荷载。
4.1 自重荷载自重荷载主要考虑连廊结构本身的重量。
根据材料的密度和结构的截面尺寸,计算出连廊结构的自重荷载。
4.2 活载荷载活载荷载主要考虑连廊上可能承载的人员和设备等活动荷载。
有限单元法课程设计有限单元法是基于连续介质力学基础上发展起来的,目前使用最广泛的数值计算方法。
有限单元法解决问题的前提是各单元相邻边界的位移协调。
有限单元法将连续的求解域离散为一组有限个单元组成的组合体,由细分单元去逼近求解域,由于单元的不同连接方式和形式各异的单元形状,因此可以适应几何形状复杂的求解区域;第二,利用每一个单元内的近似函数(形函数)来表示全求解域上待求的未知场函数,把一个连续的无限自由度问题变成离散的有限自由度问题,只要求出单元结点的物理量,就可以确定单元组合体上的其他未知场函数,如果选择合适的形函数,随着网格密度的减小,近似解将逐步趋向精确解;第三,有限单元法计算得到的总体刚度矩阵为稀疏带状矩阵,这样借助于电子计算机存储和计算的效率大大提高,便于处理大规模问题。
从上述有限单元法的特性可知,其计算原理简单,但由于单元连接方式和单元形状的多元化,以及近似函数的选择合适与否,使得有限元法在针对具体问题求解时比较烦琐,正是基于这样的应用背景,本论文提出了一种更简单实用的单元模型—平面等效桁架单元模型。
最后,编制有限元分析程序,将这种桁架单元模型运用于钢筋混凝土结构中,模型中混凝土采用等效桁架单元,钢筋采用一维杆单元,利用混凝土等效的应力应变关系对各种构件进行弹塑性分析,并试探性的提出了单元破坏准则。
用本文方法和商用有限元分析软件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,不计各杆的自重,求在荷载作用下,各杆的轴力。
作业四编制有限元程序求解结构内力位移第一部分程序设计过程和子程序的说明本例为空间桁架结构有限元分析程序。
设计思路为:自然离散桁架结构,确定各节点自由度;为单元和节点编号,输入支撑信息、荷载信息、截面特性;运行程序求得刚度矩阵,继而求得节点位移和杆件内力。
程序所需各子程序已给出,在此只需编制主程序然后调用子程序求解即可。
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. 桁架结构的原理和特点:了解桁架结构的基本原理和特点,包括节点、杆件、受力分析等内容,为后续的编程打下基础。
2. 编程软件的选择:选择合适的编程软件,例如Rhino、Grasshopper等,以及相应的插件或脚本工具,为编程工作提供支持。
3. 编程语言的学习:学习一种适合新代桁架系统编程的语言,例如Python、C#等,掌握其基本语法和常用函数,为编程工作做好准备。
三、新代桁架系统编程的实践步骤1. 数据准备:收集和整理与新代桁架系统相关的数据,包括结构参数、材料属性、受力要求等,为后续的编程工作提供输入。
2. 桁架生成:通过编程软件,根据输入的数据和设计要求,生成新代桁架系统的初始结构模型,并进行初步的受力分析。
3. 结构优化:通过编程工具,对新代桁架系统进行优化设计,以提高结构的性能和效率,例如减少材料使用量、优化节点连接等。
4. 受力分析:利用编程软件进行新代桁架系统的受力分析,验证结构的稳定性和强度,确保其满足设计要求。
5. 结果展示:通过编程软件将优化后的新代桁架系统进行可视化展示,以便设计师和工程师们更好地理解和评估结构的性能。
四、新代桁架系统编程的应用案例1. 建筑结构设计:通过新代桁架系统的编程,实现建筑结构的设计和优化,提高结构的稳定性和安全性。
2. 建筑外观设计:利用新代桁架系统的编程,实现建筑外观的设计和变形效果,为建筑增添独特的艺术魅力。
'========================================'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、变量及数组说明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)所示。