空间结构内力位移计算
第6章结构位移计算
§6-1 概述
计算构造位移旳目旳 （1）为了校核构造旳刚度。 （2）构造旳施工中，也需要构造旳位移。
图示构造进行悬臂拼装时，因为自重及吊车等荷载作用，产生位移
fA。必须先计算fA，以便采用相应措施，确保施工安全和拼装就位。
化产生旳K点旳竖向位移△Kt。α为
材料旳线膨胀系数。
ΔKt FNdut Mdt FS tds
杆轴线处旳温度变化为
t
h2 h
t1
h1 h
t2
dut tds
杆件截面对称于形心轴 t t1 t2 2
dt
t2ds t1ds
h
Δtds
h
对于杆件构造温度变化不引起剪切变形，γt=0。
t t2 t1
§6-6 静定构造温度变化时旳位移计算
将温度变化引起旳微段变形代入位移计算公式可得
ΔKt
FNtds
M
tds
h
t
FNds
t
Mds h
若各杆为等截面杆
ΔKt
tAFN
t
h
AM
AFN FN图的面积，AM M图的面积
符号旳拟定：温度变化以升温为正，轴力以拉力为正； 弯矩M以使t2边受拉为正。
对于桁架
虚拟状态如图b所示。由材料力学
dP
M Pds EI
duP
FNPds EA
Pds
kFSPds GA
k—剪切变形旳 改正系数
§6-4 静定构造在荷载作用下旳位移计算
平面杆件构造在荷载作用下旳位移计算公式为：
ΔKP
MM Pds EI
框架结构的内力和位移计算(精)
假定： （1）平面结构假定； （2）忽略柱的轴向变形； （3）D值法考虑了结点转角， 假定同层结点转角相等
2019/3/19
27
D 值法
计算方法 1、D值——修正抗侧刚度的计算 水平荷载作用下，框架不仅有侧移， 且各结点有转角，设杆端有相对位 移 ，转角 、 ，转角 1 2 位移方程为：
2019/3/19
22
反弯点法
2、剪力的计算 根据假定1：
V1 j d1 j j
Vij d ij j
Vij , d ij
——第j层第I根柱的剪力及其抗侧刚度
第j层总剪力
V pj
Vpj V1 j V2 j Vmj
2019/3/19 23
反弯点法
V1 j
第j层各柱剪力为
M ( z) N B
M(z)——上部水平荷载对坐标Z力矩总和 B——两边柱轴线间的距离
N
2019/3/19 44
柱轴向变形产生的侧移
N j
任意水平荷载下柱轴向变形产生的第j层处侧移 把框架连续化，根据单位荷载法：
2 ( NN / EA)dz
N j 0
Hj
N ( H j z) / B
框架结构的内力和位移计算荷载和设计要求51计算简图计算简图计算简图计算简图计算简图52竖向荷载作用下的近似计算方法分层法分层法分层法分层法力学知识回顾分层法计算过程构件弯矩图53水平荷载作用下内力近似计算方法反弯点法反弯点法弯点法反弯点法反弯点法反弯点法反弯点法反弯点法54水平荷载作用下内力近似计算方法d55水平荷载作用下侧移的近似计算梁柱刚度比k中柱
2019/3/19
9
计算简图
二、结构构件的截面抗弯刚度 考虑楼板的影响，框架梁的截面抗弯刚度应适当提高 现浇钢筋混凝土楼盖： 中框架：I＝2I0 边框架：I＝1.5I0 装配整体式钢筋混凝土楼盖： 截面形式选取： 框架梁跨中截面： 中框架：I＝1.5 I0 T型截面 边框架：I＝1.2 I0 框架梁支座截面： 装配式钢筋混凝土楼盖： 矩形截面 中框架：I＝I0 边框架：I＝I0 注：I0为矩形截面框架梁的截面惯性矩
框架结构的内力与位移计算
框架结构的内力与位移计算4．1 概述框架结构是目前多、高层建筑中常采用的结构形式之一。
框架在结构力学中称为刚架，结构力学中已经比较详细地介绍了超静定刚架(框架)内力和位移的计算方法，比较常用的手算方法有全框架力矩分配法、无剪力分配法和迭代法等，均为精确算法。
但在实用中大多已被更精确、更省人力的计算机分析方法(矩阵位移法)所代替。
不过，其中有些手算近似计算方法由于其计算简单、易于掌握，又能反映刚架受力和变形的基本特点，目前在实际工程中应用还很多，特别是在初步设计时的估算，手算的近似方法仍为设计人员所常用。
多、高层建筑结构在进行内力与位移计算中，为使计算简化，必须作出一些假定，以下将讨论一些结构计算中的基本假定：(1)弹性工作状态假定：结构在荷载作用下的整体工作按弹性工作状态考虑，内力和位移按弹性方法计算。
但对于框架梁及连梁等构件，可考虑局部塑性变形，内力重分布。
(2)平面结构假定：任何结构都是一个空间结构，实际风荷载及地震作用方向是随意的、不定的。
为简化计算，对规则的框架、框架—剪力墙、剪力墙结构体系及框筒结构，可将结构沿两个正交主轴方向划分为若干平面抗侧力结构—若干榀框架、若干片墙，以承受该框架、墙平面方向的水平力(风荷载及水平地震作用)，框架、墙不承受垂直于其平面方向的水平力。
(3)刚性楼面假定：各平面|考试大|抗侧力结构之间通过楼板相互联系并协同工作。
一般情况下，可认为楼板在自身平面内刚度无限大，而楼板平面外刚度很小，可以不考虑。
为保证楼面在平面内刚度，在设计中应采取相应的构造措施。
但当楼面有大开孔、楼面上有较长的外伸段、底层大空间剪力墙结构的转换层楼面以及楼面的整体性较差时，宜对采用刚性楼面假定的计算结果进行调整或在计算中考虑楼面的平面内刚度。
在上述假定下，内力分析时要解决两个问题：一个是按各片抗侧力结构的相对刚度大小，分配水平荷载至各片抗侧力结构；另一个是计算每片抗侧力结构在所分到的水平荷载作用下的内力及位移。
高层结构内力与位移计算的一般方法
a
2 1 2
；
b
1 1 2
6EJ GAL2 0
；
c
3 L0 1 2
；
e
6 1 2 L2 0
；
式中：
—截面剪应力分布不均匀系数； —中段截面的弹性模量、剪切模量、惯性矩和面积；
A E、 G 、J 、
L0 —分别为杆端刚域和中段的长度。 d j 、d k 、
j K jj Ks K kj
k K jk j K kk k
其中每个子矩阵为3×3，可叠加到总刚度矩阵 K f 中去，
叠加完毕、形成总刚度矩阵后，对于位移约束点（如底层固定点）的某一个 的位移为零时，可令相应于该位移分量的主对角元素为非常大的数（例 如 1016、 1017等），则求解时，这些位移分量自动等于零。 令荷载向量只在某一层（第i层）施加单位水平力 Pxi 1 ，其余荷载 值均为零，由下式可以求出相应的位移
GJ p , NS hNS
（1-24）
M K
（1-25）
1.1.4 在水平荷载作用下高层建筑结构空间协调工作分析的基本公式
y
us
P
( Ns ) xs
Lxs
p ys
pxs
s
vs
O′
O
lb2 hz 2 0.25hl
对于柱
lz1 hL1 0.25hz lz 2 hL 2 0.25hz
如果梁太高或太细,可按 hL hz的大小适当调整刚域长度。刚域长度当算得
为负值时，应取为零。
在进行高层建筑结构分析时，假定各片抗侧力结构在自身平面内的刚度
框架结构内力与位移计算
《高层建筑结构与抗震》辅导材料四框架结构内力与位移计算学习目标1、熟悉框架结构在竖向荷载和水平荷载作用下的弯矩图形、剪力图形和轴力图形；2、熟悉框架结构内力与位移计算的简化假定及计算简图的确定；3、掌握竖向荷载作用下框架内力的计算方法——分层法；4、掌握水平荷载作用下框架内力的计算方法——反弯点法和D值法，掌握框架结构的侧移计算方法。
学习重点1、竖向荷载作用下框架结构的内力计算；2、水平荷载作用下框架结构的内力及侧移计算。
框架在结构力学中称为刚架，刚架的内力和位移计算方法很多，可分为精确算法和近似算法。
精确法是采用较少的计算假定，较为接近实际情况地考虑建筑结构的内力、位移和外荷载的关系，一般需建立大型的代数方程组，并用电子计算机求解；近似算法对建筑结构引入较多的假定，进行简化计算。
由于近似计算简单、易于掌握，又能反映刚架受力和变形的基本特点，因此近似的计算方法仍为工程师们所常用。
本章内容主要介绍框架结构在荷载作用下内力与位移的近似计算方法。
其中分层法用于框架结构在竖向荷载作用下的内力计算，反弯点法和D值法用于框架结构在水平荷载作用下的内力计算。
既然是近似计算，就需要熟悉框架结构的计算简图和各种计算方法的简化假定。
一、框架结构计算简图的确定一般情况下，框架结构是一个空间受力体系，可以按照第四章所述的平面结构假定的简化原则，忽略结构纵向和横向之间的空间联系，忽略各构件的抗扭作用，将框架结构简化为沿横方向和纵方向的平面框架，承受竖向荷载和水平荷载，进行内力和位移计算。
结构设计时一般取中间有代表性的一榀横向框架进行分析，若作用于纵向框架上的荷载各不相同，则必要时应分别进行计算。
框架结构的节点一般总是三向受力的，但当按平面框架进行结构分析时，则节点也相应地简化。
在常见的现浇钢筋混凝土结构中，梁和柱内的纵向受力钢筋都将穿过节点或锚入节点区，这时节点应简化为刚接节点；对于现浇钢筋混凝土柱与基础的连接形式，一般也设计成固定支座，即为刚性连接。
剪力墙结构内力与位移计算3(多肢墙)
案例二：大型商业中心多肢剪力墙设计
大型商业中心多肢剪力墙设计需 要充分考虑人流、物流等动态载 荷的影响，以及商业设施的特殊
要求。
设计时应注重结构的抗震性能和 抗侧刚度，以确保在地震等突发 事件中能够保持结构的稳定性和
安全性。
设计中还需考虑商业设施的功能 需求，如商铺、餐饮、娱乐等设 施的空间布局和结构要求，以满
剪力分布
多肢剪力墙的剪力分布与墙肢的长度、 连梁的刚度以及剪力墙的整体性有关。 在墙肢和连梁的交接处，剪力通常较 大。
剪力计算公式
剪力的计算公式基于牛顿第二定律和 平衡条件，考虑剪力墙的截面特性和 荷载情况来确定。
剪力墙的轴力计算
轴力分布
多肢剪力墙的轴力分布与墙肢的长度、截面尺寸以及荷载情况有关。在墙肢和连梁的交接处，轴力可能发生突变。
足商业运营的需求。
案例三：大型桥梁多肢剪力墙设计
大型桥梁多肢剪力墙设计需要充分考虑 风、水流等自然因素的影响，以及桥梁 的跨度、荷载和施工条件等特殊要求。
设计时应采用先进的计算和分析方法， 如数值模拟或风洞试验等，对剪力墙的 受力状态进行详细分析，确保结构的稳
定性和安全性。
设计中还需考虑施工的可操作性、经济 性以及耐久性等因素，以达到最佳的设 计效果。同时，还需考虑桥梁的美观性 和环保要求，以满足社会和经济发展的
在设计多肢剪力墙时，需要充分考虑剪切传递原理，合理配置水平钢筋的位置和 数量，以确保剪力的有效传递和整体结构的稳定性。
弯曲传递原理
弯曲传递原理是指多肢剪力墙中，墙肢之间的弯矩通过竖向 钢筋传递，形成整体共同抵抗弯矩。
在设计多肢剪力墙时，需要充分考虑弯曲传递原理，合理配 置竖向钢筋的位置和数量，以确保弯矩的有效传递和整体结 构的稳定性。
结构力学大作业—多层多跨框架结构内力计算
结构力学大作业—多层多跨框架结构内力计算多层多跨框架结构是一种常见的建筑结构形式，具有良好的承载能力和空间利用率。
在设计和施工过程中，需要对该结构的内力进行计算，以确保结构的安全可靠性。
本文将介绍多层多跨框架结构内力计算的方法和步骤。
首先，需要进行结构的荷载分析。
根据设计要求和使用功能，确定结构所受的重力荷载、风荷载和地震荷载等。
根据荷载的作用位置和方向，进行荷载分解和转换，将荷载分配给各个构件。
其次，进行结构模型的建立。
采用有限元分析软件（如ANSYS、ABAQUS等）或者手工计算的方法，将结构划分为有限个构件，建立相应的杆件模型。
注意考虑构件的截面形状、材料性质和连接方式等。
然后，进行结构的静力分析。
根据结构的平衡条件和变形规律，建立结构的刚性方程组。
通过求解刚性方程组，得到结构的各个节点的位移和构件内力。
需要注意的是，在多层多跨框架结构中，由于存在多个层和多个跨度，结构的刚性方程组会比较复杂，需要采用适当的计算方法（如直接刚度法、位移法等）进行求解。
最后，进行结构的动力分析。
根据结构的动力特性（如固有周期、振型等），确定结构所受的地震荷载或者其他动力荷载，并进行结构的动力响应分析。
通过求解结构的动力方程组，得到结构的振动位移和构件的动力响应。
在进行内力计算时，需要考虑结构的整体受力平衡和构件的局部受力平衡。
对于整体受力平衡，需要保证结构的静力平衡，即各个节点的受力之和为零。
对于局部受力平衡，需要考虑构件的内力传递和受力集中等因素，以确保构件的强度和稳定性。
在实际计算过程中，需要根据具体的结构形式和工程要求，选择适当的计算方法和假设条件。
特别是在多层多跨框架结构中，由于受力情况的复杂性，需要进行合理的简化和近似，以提高计算效率和准确性。
综上所述，多层多跨框架结构内力计算是结构力学领域的重要内容之一、通过对结构的荷载分析、模型建立、静力分析和动力分析等步骤的综合运用，可以得到结构的内力分布情况，为结构的设计和施工提供科学依据。
空间结构中的力学行为分析
空间结构中的力学行为分析在建筑和工程领域中，空间结构是指由梁、柱、墙等构件组成的三维结构体系。
这些构件在承受荷载时会发生力学行为，因此对空间结构的力学行为进行分析是非常重要的。
一、力学行为的分类空间结构的力学行为可以分为静力学行为和动力学行为两类。
静力学行为是指在静力平衡条件下，结构在荷载作用下的变形和力的分布情况。
静力学行为的分析是建筑和工程设计的基础，通过计算和模拟，可以确定结构的稳定性和安全性。
动力学行为是指结构在受到外界激励时的响应情况。
这些激励可以是地震、风力、交通振动等。
动力学行为的分析可以帮助设计师评估结构的抗震性能和振动特性，从而进行合理的结构设计和改进。
二、力学行为的分析方法1. 静力学分析方法静力学分析是通过建立结构的数学模型，应用静力学原理和力学方程进行计算和分析。
常用的静力学分析方法包括刚度法、弹性理论、有限元法等。
刚度法是一种基于刚度的分析方法，通过计算结构的刚度矩阵和荷载向量，可以求解结构的位移和内力。
弹性理论基于弹性材料的假设，将结构看作由弹性构件组成的整体，通过应力-应变关系来分析结构的响应。
有限元法则是将结构离散化为有限个单元，通过求解节点的位移和应力，来模拟整个结构的力学行为。
2. 动力学分析方法动力学分析是通过建立结构的动力学模型，考虑结构的质量、刚度和阻尼等因素，来研究结构在外界激励下的响应。
常用的动力学分析方法包括模态分析、响应谱分析和时程分析等。
模态分析是一种基于结构的固有振动模态进行分析的方法。
通过求解结构的固有频率、振型和阻尼比等参数，可以评估结构的振动特性和抗震性能。
响应谱分析是一种基于结构的响应谱进行分析的方法，通过结合地震动的频谱和结构的频率响应，来评估结构在地震作用下的最大响应。
时程分析是一种基于结构的动力学方程进行数值求解的方法，通过模拟结构在时间上的响应，来研究结构的振动行为。
三、力学行为的影响因素空间结构的力学行为受到多种因素的影响，包括结构的几何形状、材料的性质、荷载的大小和方向等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
作业四编制有限元程序求解结构内力位移第一部分程序设计过程和子程序的说明本例为空间桁架结构有限元分析程序。
设计思路为：自然离散桁架结构，确定各节点自由度；为单元和节点编号，输入支撑信息、荷载信息、截面特性；运行程序求得刚度矩阵，继而求得节点位移和杆件内力。
程序所需各子程序已给出，在此只需编制主程序然后调用子程序求解即可。
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）程序运行完毕，查看输出数据文件结果。
2、参照附录部分各个变量的含义和说明，填写输入文件（1）输入NP,NE,NR,NM,NCF（2）输入X(M),Y(M),Z(M),M=1,NP（3）输入MM(M),ME(1,M),ME(2,M),M=1,NE（4）输入RR(M),M=1,NR（5）输入PF(1,M),PF(2,M),PF(3,M),PF(4,M),M=1,NCF（6）输入AE(1,1),AE(2,1)7第三部分空间桁架有限元程序1、源程序PROGRAM MAINOPEN(1,FILE="IN.TXT",STATUS='OLD')OPEN(2,FILE="OUT1.TXT",STATUS='UNKNOWN')READ(1,*)NP,NE,NR,NM,NCFCALL LAST(NP,NE,NR,NM,NCF)ENDSUBROUTINE LAST(NP,NE,NR,NM,NCF)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4(I-N)DIMENSION X(NP),Y(NP),Z(NP),MM(NE),ME(2,NE),&RR(2,NR),PF(4,NCF),AE(2,NM),SG(NE),SM(NE),&IT(3,NP),DIST(NP*3),LMT(6,NE),MAXA(1000),V(1000),CKK(1000),& FF(NP*3),PP(NP*3)READ(1,*)(X(M),Y(M),Z(M),M=1,NP)READ(1,*)(MM(I),ME(1,I),ME(2,I),I=1,NE)READ(1,*)(RR(1,M),RR(2,M),M=1,NR)READ(1,*)(PF(1,M),PF(2,M),PF(3,M),PF(4,M),M=1,NCF)READ(1,*)(AE(1,M),AE(2,M),M=1,NM)WRITE(2,'(20X,A)')'///////////////////////////////////////////'WRITE(2,'(20X,A)')'/ program for analysis of space trusses /'WRITE(2,'(20X,A)')'/ school of civil engineering CSU /'WRITE(2,'(20X,A)')'/ 2012.06.23 designed by HuLihua /'8WRITE(2,'(20X,A/)')'///////////////////////////////////////////'WRITE(2,'(20X,/A/)')'the control data:'WRITE(2,10)'num joints','num members','num supports','num type','num p' 10 FORMAT(1X,A,2X,A,2X,A,4X,A,4X,A)WRITE(2,15)NP,NE,NR,NM,NCF15 FORMAT(3X,I4,8X,I4,9X,I4,8X,I5,8X,I5)WRITE(2,'(/A/)')'the information of joints:'WRITE(2,70)'num','x','y','z'70 FORMAT(1X,A,9X,A,9X,A,10X,A)WRITE(2,80)(I,X(I),Y(I),Z(I),I=1,NP)80 FORMAT(1X,I4,3X,E8.2,3X,E8.2,3X,E8.2)WRITE(2,'(/A/)')'the information of members:'WRITE(2,25)'member','start','end'25 FORMAT(1X,A,2X,A,3X,A)WRITE(2,30)(MM(I),ME(1,I),ME(2,I),I=1,NE)!WRITE(*,*)(ME(1,I),I=1,NE)30 FORMAT(1X,I4,3X,I4,3X,I4)WRITE(2,'(/A/)')'the information of supports:'WRITE(2,35)'joint','s'35 FORMAT(1X,A,2X,A)WRITE(2,400)(RR(1,M),RR(2,M),M=1,NR)400 FORMAT(1X,F4.0,2X,F5.3)WRITE(2,'(/A/)')'the information of load:'WRITE(2,900)'member','x','y','z'9900 FORMAT(1X,A,2X,A,3X,A,8X,A)WRITE(2,600)(PF(1,M),PF(2,M),PF(3,M),PF(4,M),M=1,NCF)600 FORMAT(1X,F8.0,3X,F8.3,3X,F8.3,3X,F8.3)WRITE(2,'(/A/)')'the information of sections:'WRITE(2,700)'a','e'700 FORMAT(1X,A,2X,A)WRITE(2,800)AE(1,1),AE(2,1)800 FORMAT(1X,1PE8.2,3X,F8.4)CALL FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT)CALL FMAXA(NN1,NE,LMT,MAXA,NWK,NP)!WRITE(*,*)(LMT(1,I),LMT(2,I),LMT(3,I),LMT(4,I),LMT(5,I),LMT(6,I) ,I=1,26)!WRITE(*,*)(MAXA(I),I=1,NN1)CALL CONKB(NP,NE,NWK,ME,X,Y,Z,AE,LMT,MAXA,CKK,NN1)!WRITE(*,*)(TAK(I,1),TAK(I,2),TAK(I,3),TAK(I,4),TAK(I,5),TAK(I,6), I=1,6)!WRITE(*,*)NWK!WRITE(*,*)(CKK(I),I=1,NWK)!CALL P1(NCF,P,PF,NP)!CALL SUPPORT(NP,NR,JS,IDS,DJS,TK,P)ISH=10CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1)CALL REBACK(CKK,V,MAXA,NN,NWK,NN1,NP,NCF,PF,IT)WRITE(2,'(20X,///A/)')'the result data:'WRITE(2,'(/1X,A/)')'joint displacements'10WRITE(2,'(1X,A,6X,A,15X,A,15X,A/)')'joint','U','V','W'!WRITE(2,40)(I,P(3*I-2),P(3*I-1),P(3*I),I=1,NP)!40 FORMAT(I3,4X,E11.4,4X,E11.4,4X,E11.4)!DO IE=1,NECALL DISPLS(NP,NE,NN,IT,V,DIST,AE,X,Y,Z,PP,FF,SG,SM,ME) WRITE(2,40)(I,DIST(3*I-2),DIST(3*I-1),DIST(3*I),I=1,NP)40 FORMAT(I3,4X,E11.4,4X,E11.4,4X,E11.4)!END DOWRITE (2,'(20X,/A/)')'axial force'WRITE(2,'(/1X,A,5X,A/)')'member','N'WRITE(2,50)(M,SG(M),M=1,NE)50 FORMAT(I3,4X,E10.2)WRITE (2,'(20X,/A/)')'stress'WRITE(2,'(/1X,A,5X,A/)')'member','σ'WRITE(2,60)(M,SM(M),M=1,NE)60 FORMAT(I3,4X,E10.2)RETURNEND!矩阵转置SUBROUTINE MAT(M,N,A,B)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4(I-N)DIMENSION A(M,N),B(N,M)DO I=1,M11DO J=1,NB(J,I) = A(I,J)END DOEND DORETURNEND!形成单元局部坐标系中的刚度矩阵SUBROUTINE FKE(IE,NP,NE,X,Y,Z,ME,AE,AKE)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4 (I-N)DIMENSION AKE(2,2),X(NP),Y(NP),Z(NP),ME(2,NE),AE(2,1) !WRITE(*,*)ME(1,IE)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)WRITE(*,*)BLE0=AE(1,1)A0=AE(2,1)C=E0*A0/BL!WRITE(*,*)CAKE(1,1)=CAKE(1,2)=-C12AKE(2,1)=-CAKE(2,2)=CRETURNEND!局部坐标与总体坐标转换矩阵SUBROUTINE FT(IE,NP,NE,X,Y,Z,ME,T)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4 (I-N)DIMENSION X(NP),Y(NP),Z(NP),ME(2,NE),T(2,6) !WRITE(*,*)NP!WRITE(*,*)ME(1,IE)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)!WRITE(*,*)X1,X2BL=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)=CZ!WRITE(*,*)(T(1,I),T(2,I),I=1,6)13RETURNEND!形成IT和LMT数组SUBROUTINE FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4(I-N)!This program forms the joint&element numbering matrix IT&LMTDIMENSION IT(3,NP),LMT(6,NE),ME(2,NE),RR(2,NR)NN=0;NN1=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=CC=C-NCDO J=1,3C=C*10.0L=C+0.1C=C-LIF(L.EQ.0)THENN=N+114IT(J,I)=NELSEIT(J,I)=0ENDIFENDDOENDDONN=NNN1=NN+1DO IE=1,NEDO I=1,2NI=ME(I,IE)DO J=1,3LMT((I-1)*3+J,IE)=IT(J,NI) ENDDOENDDOENDDORETURNEND!矩阵的乘法SUBROUTINE MUL(A,B,M,N,L,AB) IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4 (I-N) DIMENSION A(M,N),B(N,L),AB(M,L) DO I=1,M15DO J=1,LAB(I,J)=0DO K=1,NAB(I,J)=AB(I,J)+A(I,K)*B(K,J)ENDDOENDDOENDDOEND!形成MAXA数组SUBROUTINE FMAXA(NN1,NE,LMT,MAXA,NWK,NP) IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4 (I-N)! This program forms the MDE address matrix MAXA of [K] DIMENSION MAXA(3*NP),LMT(6,NE)MAXA=0;NWK=0MAXA(1)=1DO I=2,NN1IP=I-1IG=IPDO IE=1,NEDO J=1,6IF(LMT(J,IE).EQ.IP) THENDO K=1,6IF(LMT(K,IE).GT.0.AND.LMT(K,IE).LE.IG) IG=LMT(K,IE)16ENDDOEND IFENDDOENDDOMAXA(I)= MAXA(I-1)+IP-IG+1ENDDONWK=MAXA(NN1)-1RETURNEND!形成一维存储数组总体刚度矩阵CKKSUBROUTINECONKB(NP,NE,NWK,ME,X,Y,Z,AE,LMT,MAXA,CKK,NN1)IMPLICIT REAL*8 (A-H,O-Z)IMPLICIT INTEGER*4 (I-N)DIMENSION CKK(NWK),X(NP),Y(NP),Z(NP),AE(2,1),&LMT(6,NE),ME(2,NE),MAXA(NN1),&AKE(2,2),T(2,6),TT(6,2),TAK(6,6),AK(6,2)!WRITE(*,*)(LMT(1,I),LMT(2,I),LMT(3,I),LMT(4,I),LMT(5,I),L MT(6,I),I=1,26)CKK=0!WRITE(*,*)(ME(1,IE),IE=1,NE)DO IE=1,NE!WRITE(*,*)NETAK=0K=IE17CALL FKE(IE,NP,NE,X,Y,Z,ME,AE,AKE)CALL FT(IE,NP,NE,X,Y,Z,ME,T)!WRITE(*,*)(T(1,I),T(2,I),I=1,6)CALL MAT(2,6,T,TT)CALL MUL(TT,AKE,6,2,2,AK)CALL MUL(AK,T,6,2,6,TAK)!CALL FLMT(NP,NE,NN,NN1,NR,RR,ME,IT,LMT) !CALL FMAXA(NN1,NE,LMT,MAXA,NWK,NP)!WRITE(*,*)LMT(1,IE)DO I=1,6DO J=1,6NI=LMT(I,IE)!WRITE(*,*)NINJ=LMT(J,IE)IF((NJ-NI).GE.0.AND.NI*NJ.GT.0)THENK=MAXA(NJ)+NJ-NICKK(K)=CKK(K)+TAK(I,J)ENDIFENDDOENDDOENDDORETURNEND!解线性方程组得位移值V18!SUBROUTINELDLTREBACK(PF,CKK,V,MAXA,IT,NN,NWK,NN1,ISH,IOUT,NCF,N P)!IMPLICIT REAL*8(A-H,O-Z)!IMPLICIT INTEGER*4(I-N)!DIMENSIONCKK(NWK),V(NN,1),MAXA(NN1),PF(4,NCF),IT(3,NP)!CALL LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1)!CALL REBACK(CKK,V,MAXA,NN,NWK,NN1,NP,NCF,PF,IT)!END!将一维存储的结构刚度矩阵进行LDLT分解SUBROUTINE LDLT(CKK,MAXA,NN,ISH,IOUT,NWK,NN1)IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4 (I-N)DIMENSION CKK(NWK),MAXA(NN1)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=KU19DO 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+CKK(KI+L)*CKK(KLT+L)CKK(KLT)=CKK(KLT)-C260 K=K+1240 K=NB=0.0DO 300 KK=KL,KUK=K-1KI=MAXA(K)C=CKK(KK)/CKK(KI)IF(ABS(C).LT.1.0E+07) GOTO 290WRITE(IOUT,2010) N,CSTOP290 B=B+C*CKK(KK)300 CKK(KK)=CCKK(KN)=CKK(KN)-B20304 IF(CKK(KN)) 310,310,200310 IF(ISH.EQ.0) GOTO 320IF(CKK(KN).EQ.0.0) CKK(KN)=-1.0E-16GOTO 200320 WRITE(IOUT,2000) N,CKK(KN)STOP200 CONTINUERETURN2000 FORMAT(//' Stop-stiffness matrix not positive& + definite',//,'nonpositive pivot for equation',&I4,//,' pivot =',E20.10)2010 FORMAT(//,' Stop-sturm sequence check failed& + because of multiplier growth for column&+ number',I4,//, ' Multiplier = ',E20.8)END!回代求解子程序SUBROUTINEREBACK(CKK,V,MAXA,NN,NWK,NN1,NP,NCF,PF,IT) IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4 (I-N)DIMENSIONCKK(NWK),V(NN),MAXA(NN1),PF(4,NCF),IT(3,NP) V=0DO J=1,NCFM=PF(1,J)21L=IT(1,M)V(L)=PF(2,J)V(L+1)=PF(3,J)V(L+2)=PF(4,J)ENDDODO 400 N=1,NNKL=MAXA(N)+1KU=MAXA(N+1)-1IF(KU-KL) 400,410,410 410 K=NC=0.0DO 420 KK=KL,KUK=K-1420 C=C+CKK(KK)*V(K)V(N)=V(N)-C400 CONTINUEDO 480 N=1,NNK=MAXA(N)480 V(N)=V(N)/CKK(K)IF(NN.EQ.1)RETURNN=NNDO 500 L=2,NNKL=MAXA(N)+1KU=MAXA(N+1)-122IF(KU-KL) 500,510,510510 K=NDO 520 KK=KL,KUK=K-1520 V(K)=V(K)-CKK(KK)*V(N)500 N=N-1RETURNEND!求单元内力节点位移及约束SUBROUTINEDISPLS(NP,NE,NN,IT,V,DIST,AE,X,Y,Z,PP,FF,SG,SM,ME) IMPLICIT REAL*8(A-H,O-Z)IMPLICIT INTEGER*4(I-N)DIMENSIONIT(3,NP),DIST(NP*3),V(NN),T(2,6),X(NP),Y(NP),Z(NP),&TT(6,2),AE(2,1),ME(2,NE),UE(6),U(2),&AKE(2,2),FE1(2),FE(6),FF(NP*3),PP(NP*3),SG(NE),SM(NE) SG=0;SM=0;FF=0DO I=1,NPDO J=1,3LAB=IT(J,I)IF(LAB.EQ.0) THENDIST(3*(I-1)+J)=0.0ELSEIF(B.LE.NN) THENDIST(3*(I-1)+J)=V(LAB)23ENDIFENDDOENDDODO IE=1,NEN1=ME(1,IE);N2=ME(2,IE)DO J=1,3UE(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(IE,NP,NE,X,Y,Z,ME,AE,AKE)CALL MUL(T,UE,2,6,1,U) !U=T*UE求局部坐标系中的节点位移CALL MUL(AKE,U,2,2,1,FE1) !FE1=AKE*U求单元局部坐标系节点力CALL MAT(2,6,T,TT)CALL MUL(TT,FE1,6,2,1,FE) !FE=TT*FE1DO J=1,3FF(3*(N1-1)+J)=FF(3*(N1-1)+J)+FE(J)FF(3*(N2-1)+J)=FF(3*(N2-1)+J)+FE(3+J)ENDDOAO=AE(2,1)SG(IE)=FE1(2)SM(IE)=FE1(2)/AODO I=1,NP*32425FF(I)=FF(I)-PP(I) ENDDO ENDDO RETURN END第四部分 程序使用算例121110987654321Y ZX2019181716151413121110823765419444421222324252610kNkN10说明：A=0.002，E=2.1E8输出文件OUT1 //////////////////////////////////////////// program for analysis of space trusses / / school of civil engineering CSU / / 2012.06.23 designed by Hulihua ////////////////////////////////////////////the control data:num joints num members num supports num tyepe num p12 26 4 1 2 the information of joints:num x y z1 .00E+00 .00E+00 .00E+002 .40E+01 .00E+00 .00E+003 .80E+01 .00E+00 .00E+004 .00E+00 .40E+01 .00E+005 .40E+01 .40E+01 .00E+006 .80E+01 .40E+01 .00E+007 .00E+00 .00E+00 .40E+018 .40E+01 .00E+00 .40E+019 .80E+01 .00E+00 .40E+0110 .00E+00 .40E+01 .40E+0111 .40E+01 .40E+01 .40E+0112 .80E+01 .40E+01 .40E+01the information of members:member start end1 1 22 2 33 1 44 2 55 3 66 4 57 5 68 7 8269 8 910 7 1011 8 1112 9 1213 10 1114 11 1215 1 816 3 817 4 1118 6 1119 1 1020 3 1221 1 722 2 823 3 924 4 1025 5 1126 6 12the information of supports:joint s1. .1113. .1114. .1116. .111the information of load:member x y z2710. .000 10.000 .00012. .000 10.000 .000the information of sections:a e2.10E+08 .02000the result data:joint displacementsjoint U V W1 .0000E+00 .0000E+00 .0000E+002 .0000E+00 .0000E+00 .0000E+003 .0000E+00 .0000E+00 .0000E+004 .0000E+00 .0000E+00 .0000E+005 .0000E+00 .0000E+00 .0000E+006 .0000E+00 .0000E+00 .0000E+007 .0000E+00 .3646E-03 .0000E+008 .0000E+00 .0000E+00 .0000E+009 .0000E+00 .3646E-03 .0000E+0010 .0000E+00 .3646E-03 -.9524E-0411 .0000E+00 .0000E+00 .0000E+0012 .0000E+00 .3646E-03 -.9524E-04axial forcemember N1 .00E+002 .00E+003 .00E+00285 .00E+006 .00E+007 .00E+008 .00E+009 .00E+0010 .00E+0011 .00E+0012 .00E+0013 .00E+0014 .00E+0015 .00E+0016 .00E+0017 .00E+0018 .00E+0019 .14E+0220 .14E+0221 .00E+0022 .00E+0023 .00E+0024 -.10E+0225 .00E+0026 -.10E+02stressmember σ1 .00E+00293 .00E+004 .00E+005 .00E+006 .00E+007 .00E+008 .00E+009 .00E+0010 .00E+0011 .00E+0012 .00E+0013 .00E+0014 .00E+0015 .00E+0016 .00E+0017 .00E+0018 .00E+0019 .71E+0420 .71E+0421 .00E+0022 .00E+0023 .00E+0024 -.50E+0425 .00E+0026 -.50E+0430附录程序各种变量的说明变量含义及说明控制数据NE 杆件单元数NP 节点数NR 支撑条件数NM 截面特征类型数NCF 集中荷载数目结构描述数据ME(1,NE) 杆件始端节点号ME(2,NE) 杆件末端节点号X(NP) 节点X坐标Y(NP) 节点Y坐标Z(NP) 节点Z坐标RR(1,NP) 支撑条件所处节点号RR(2,NP) 约束条件PF(1,NP) 施加荷载的节点号PF(2,NP) X方向上施加荷载的大小（KN）PF(3,NP) Y方向上施加荷载的大小（KN）PF(4,NP) Z方向上施加荷载的大小（KN）AE(1,1) 弹性模量EIT（3，NP）节点联系数组AE(2,1) 截面尺寸AAKE 单元在局部坐标系中的刚度矩阵TAK 单元在整体坐标系中的刚度矩阵LMT 单元联系数组T 坐标转换矩阵MAXA(NN1) 整体坐标系中的刚度矩阵对角线自由度数组CKK 一维存储总体刚度矩阵SG 杆件内力SM 杆件应力DIST(3*I-2) 单元I节点X方向的位移DIST(3*I-1) 单元I节点Y方向的位移DIST(3*I) 单元I节点Z方向的位移BL 杆件长度NN 结构矩阵方程的方程总数NWK 一维存储的结构刚度矩阵元素总数31。