平面桁架程序计算原理及程序编制

1、平面桁架程序计算原理及程序编制 结构有限元程序设计基本原理 平面桁架程序的计算原理及程序编制平面桁架程序的计算原理及程序编制 平面桁架程序计算原理及程序编制 5-01 矩阵位移法矩阵位移法 5-02 矩阵位移法算例矩阵位移法算例 5-03 最小总势能原理的应用最小总势能原理的应用 5-04 矩阵位移法的求解步骤矩阵位移法的求解步骤 5-05 结构计算简图的数据结构结构计算简图的数据结构 5-06 位移未知数的确定位移未知数的确定 5-07 单根杆件的分析单根杆件的分析 5-08 结构总刚度矩阵的形成结构总刚度矩阵的形成 5-09 杆件内力的计算杆件内力的计算 5-10 能量原理和矩阵位移法能量

2、原理和矩阵位移法 平面桁架程序的计算原理平面桁架程序的计算原理 平面桁架程序计算原理及程序编制 5-1 5-1 矩阵位移法矩阵位移法 桁架是由离散杆件组成的构架结构，杆件的端点桁架是由离散杆件组成的构架结构，杆件的端点 借助于无摩擦的铰连接起来。借助于无摩擦的铰连接起来。 桁架主要靠各杆中的轴向拉力或压力来传递作用桁架主要靠各杆中的轴向拉力或压力来传递作用 在桁架节点上的荷载，杆件的任何弯曲均忽略不在桁架节点上的荷载，杆件的任何弯曲均忽略不 计。计。 用有限元分析桁架时，桁架中的每根杆件都是一用有限元分析桁架时，桁架中的每根杆件都是一 个单元，称为杆单元。它是一维单元。个单元，称为杆单元。它是

3、一维单元。 不同分类：不同分类： （1 1）平面桁架）平面桁架/ /空间桁架空间桁架 （2 2）静定桁架）静定桁架/ /超静定桁架超静定桁架 平面桁架程序计算原理及程序编制 求解方法 力法：以力未知数，必须预先满足平衡条件，然力法：以力未知数，必须预先满足平衡条件，然 后通过连续条件求解未知力；超静定基的选取。后通过连续条件求解未知力；超静定基的选取。 位移法：以位移为未知数，各杆件的变形由相连位移法：以位移为未知数，各杆件的变形由相连 接的节点位移确定（变形协调条件），通过各个接的节点位移确定（变形协调条件），通过各个 节点的平衡方程求出未知位移，再由位移计算出节点的平衡方程求出未知位移，再

4、由位移计算出 各杆件的内力；各节点的平衡方程也可由最小总各杆件的内力；各节点的平衡方程也可由最小总 势能原理推导。势能原理推导。 以平面桁架结构分析的程序设计为例，介绍结构分以平面桁架结构分析的程序设计为例，介绍结构分 析和程序设计的方法。析和程序设计的方法。 平面桁架程序计算原理及程序编制 5-2 5-2 矩阵位移法算例矩阵位移法算例 如图所示平面桁架，各根杆的截面积如图所示平面桁架，各根杆的截面积F F相等，材料的相等，材料的 弹性模量弹性模量E E相同，在两个单位力相同，在两个单位力P=1P=1的外荷载作用下，的外荷载作用下， 用位移法计算节点位移和各杆内力。用位移法计算节点位移和各杆内

5、力。 平面桁架结构 12 34 PP 节点总数：节点总数： NW=4NW=4 可动节点数：可动节点数：NU=2NU=2 位移未知数总数：位移未知数总数：NDISP=2NDISP=2* *NU=4NU=4 1 1 2 2 u v W u v 注意：节点自由注意：节点自由 度排序和节点平度排序和节点平 衡方程相对应。衡方程相对应。 平面桁架程序计算原理及程序编制 根据虎克定律，对于任意一根杆，有：根据虎克定律，对于任意一根杆，有： E F NL L 考虑单根杆件在位移下产生的内力。考虑单根杆件在位移下产生的内力。 式中式中 为杆的伸长，为杆的伸长， 为杆的长度，为杆的长度， 为杆的内为杆的内 力，

6、力， 称为单根杆的刚度称为单根杆的刚度( (单元刚度阵单元刚度阵) )。 LL N EF L 单元局部刚度阵单元局部刚度阵 用用RDRD来表示单元刚度阵，于是来表示单元刚度阵，于是 EF RD L 胡克定律可表达如下胡克定律可表达如下 NRDL 平面桁架程序计算原理及程序编制 由单根杆件的变形几何关系可得由单根杆件的变形几何关系可得 cossin () cos() sin BABA Luv uuvv 角是杆件轴线与角是杆件轴线与 方向的夹角，由方向的夹角，由 正向逆时正向逆时 针向转至杆轴的角度为正。针向转至杆轴的角度为正。 xx A A u v B B u v A B x y 进一步有内力和

7、节点位移之间的关系进一步有内力和节点位移之间的关系 () cos() sin BABA EF Nuuvv L 平面桁架程序计算原理及程序编制 对于每根杆件（以两端节点编号对于每根杆件（以两端节点编号A A和和B B定出定出 角）角） 应用上述公式，有应用上述公式，有 1 1 v NEF a 12 2 uu NEF a 2 3 v NEF a 11 4 2 uv NEF a 21 5 2 vu NEF a 应用上述各杆内力和位移关系后，便可建立以位应用上述各杆内力和位移关系后，便可建立以位 移为未知数的节点平衡方程。移为未知数的节点平衡方程。 平面桁架结构 平面桁架程序计算原理及程序编制 节点平

8、衡方程为节点平衡方程为 0 24 cos450NNP 0 14 cos450NN 0 25 cos450NNP 0 35 cos450NN 注意：节点自由度排序和节点平衡方程相对应。注意：节点自由度排序和节点平衡方程相对应。 平面桁架程序计算原理及程序编制 将各杆的内力用位移表示的方程代入上式，有将各杆的内力用位移表示的方程代入上式，有 112 11 10 2 22 2 EF uvuP a 11 11 100 2222 EF uvP a 122 11 01 2 22 2 EF uuvP a 22 11 001 2 22 2 EF uvP a 未知未知 数为数为 节点节点 位移位移 平面桁架程序

9、计算原理及程序编制 写成矩阵形式写成矩阵形式 RWP，即正则方程。，即正则方程。 其中其中 11121314 21222324 31323334 41424344 11 110 2 22 2 11 100 2 22 2 / 11 101 2 22 2 11 001 2 22 2 rrrr rrrr REF a rrrr rrrr 结构刚度矩阵结构刚度矩阵 而外力向量为：而外力向量为： 1 2 3 4 0 0 PP P P PP P 平面桁架程序计算原理及程序编制 采用高斯消元法求得位移为采用高斯消元法求得位移为 1 3.824 /ua EF 1 1.000 /va EF 21 uu 21 vv

10、 利用每根杆的内力利用每根杆的内力- -位移关系计算杆内力位移关系计算杆内力 1 1N 2 0N 3 1N 4 1.412N 5 1.412N 平面桁架程序计算原理及程序编制 5-3 5-3 最小总势能原理的应用最小总势能原理的应用 总势能由两部分组成总势能由两部分组成 结构的弹性应变能结构的弹性应变能 外力由于结构变位所产生的势能外力由于结构变位所产生的势能 U V n最小总势能原理与位移法都是以位移为未知最小总势能原理与位移法都是以位移为未知 数使变形状态预先满足连续条件。现对上述数使变形状态预先满足连续条件。现对上述 例题采用最小总势能原理进行求解。例题采用最小总势能原理进行求解。 平面

11、桁架程序计算原理及程序编制 对于整个桁架应变能是所有杆件应变能的叠对于整个桁架应变能是所有杆件应变能的叠 加，即加，即 5 2 1 1 () () / 2 kkk k UEFLL 由于桁架结构的杆件只能承受拉压力，所以由于桁架结构的杆件只能承受拉压力，所以 单根杆件的应变能为单根杆件的应变能为 2 () / 2 EF ULL 平面桁架程序计算原理及程序编制 把上述值代入应变能表达式，得到把上述值代入应变能表达式，得到 22222 11221222 ()() /2 2() /2 2 2 EF Uvuuvuvvu a 由上式可见，公式中只有位移的二次项，也就是说由上式可见，公式中只有位移的二次项，

12、也就是说 是位移的一个二次齐次函数，或者说是一个位移的是位移的一个二次齐次函数，或者说是一个位移的 二次型。（位移的正定二次型，应变能总是正的）二次型。（位移的正定二次型，应变能总是正的） 对于该平面桁架，有对于该平面桁架，有 杆号杆号12345 Laaasqrt(2)asqrt(2)a Lv1u1-u2v2(u1+ v1)/ sqrt(2)(v2- u2)/ sqrt(2) 平面桁架程序计算原理及程序编制 现在计算外力势能。外力现在计算外力势能。外力 产生势能的原因是当节产生势能的原因是当节 点发生位移时，外力要作功。所作功的负值便是它点发生位移时，外力要作功。所作功的负值便是它 们具有势能

13、的改变量，如果取未变形位置外力的势们具有势能的改变量，如果取未变形位置外力的势 能为零，有能为零，有 P 12 VPuPu 平面桁架结构 12 34 PP 将将 和和 相加，得到总势能。相加，得到总势能。 由于由于 是位移的一次函数，总势能就成为位是位移的一次函数，总势能就成为位 移的二次非齐次函数。移的二次非齐次函数。 UV V 平面桁架程序计算原理及程序编制 根据最小总势能原理，在所有可能的位移状态中，根据最小总势能原理，在所有可能的位移状态中， 真正发生的位移状态使总势能为最小。即函数对真正发生的位移状态使总势能为最小。即函数对 自变量的偏微商为零，即自变量的偏微商为零，即 () 0(1

14、,2,3,4) i UV i w 11 21 32 42 wu wv W wu wv 式中式中 平面桁架程序计算原理及程序编制 现对各位移变量分别取偏微商后，得现对各位移变量分别取偏微商后，得 112 11 10 2 22 2 EF uvuP a 11 11 100 2222 EF uvP a 122 11 01 2 22 2 EF uuvP a 22 11 001 2 22 2 EF uvP a 注：值得指出的是刚度矩阵中的系数只与结构本身的几何形态和约束注：值得指出的是刚度矩阵中的系数只与结构本身的几何形态和约束 条件有关，而与外载无关。条件有关，而与外载无关。 平面桁架程序计算原理及程序

15、编制 5-4 5-4 矩阵位移法的求解步骤矩阵位移法的求解步骤 u结构计算简图结构计算简图-（节点、单元编号，建立一个统一的坐标（节点、单元编号，建立一个统一的坐标 系等）系等） u分析节点位移的力学特性分析节点位移的力学特性-确定位移未知数。确定位移未知数。 u建立每根杆件两端位移和内力的关系建立每根杆件两端位移和内力的关系-（单根杆件的刚度（单根杆件的刚度 矩阵）矩阵） u根据每根杆件上的上述关系建立结构可动节点的平衡方程根据每根杆件上的上述关系建立结构可动节点的平衡方程 -（结构总刚度矩阵）（结构总刚度矩阵） u求解平衡方程，得到节点位移。求解平衡方程，得到节点位移。 u根据求得的位移，

16、利用每根杆件位移与内力关系计算各杆根据求得的位移，利用每根杆件位移与内力关系计算各杆 的内力。的内力。 平面桁架程序计算原理及程序编制 5-5 5-5 结构计算简图的数据结构结构计算简图的数据结构 完整而确切描述一个平面桁架结构的数据有三个方面：完整而确切描述一个平面桁架结构的数据有三个方面： （1 1）结构本体描述数据（）结构本体描述数据（NW, IESG, NU, X, Y, HL, HRNW, IESG, NU, X, Y, HL, HR） （2 2）性质数据（）性质数据（EFEF） （3 3）荷载数据（）荷载数据（P PX X, P, PY Y） NWNW为节点总数为节点总数 IESG

17、IESG杆件总数杆件总数 NUNU可动节点总数可动节点总数 X, Y X, Y 节点坐标节点坐标 HL, HR HL, HR 每根杆件两端节点编号每根杆件两端节点编号 EF EF 性质数据性质数据 P PX X, P, PY Y 外载荷数据 外载荷数据 平面桁架程序计算原理及程序编制 5-6 5-6 位移未知数的确定位移未知数的确定 对于平面桁架，每个节点有两个自由对于平面桁架，每个节点有两个自由 度，把第个度，把第个 节点的水平位移、垂直节点的水平位移、垂直 位移分别记为位移分别记为 ， ，这样结构共有，这样结构共有 2NU2NU个位移。个位移。 位移的方向与坐标轴相同为正，以这位移的方向与

18、坐标轴相同为正，以这 些位移作为未知数，并排列成一个列些位移作为未知数，并排列成一个列 向量，称为结构总位移向量。向量，称为结构总位移向量。 i i u i v 11 12 23 2 i i nu nuNDISP uw vw uw v W u v u vw 第第i个节点的位移在总位移向量中占个节点的位移在总位移向量中占 i0+1和和i0+2位置，而位置，而i0=2*i -2。 平面桁架程序计算原理及程序编制 5-7 5-7 单根杆件的分析单根杆件的分析 根据虎克定律，对于任意一根杆，有：根据虎克定律，对于任意一根杆，有： EF NL L 式中式中 为杆的伸长，为杆的伸长， 为杆的长度，为杆的长

19、度， 为杆的为杆的 内力，内力， 称为单根杆的刚度（单元刚度阵）。称为单根杆的刚度（单元刚度阵）。 LLN /EFL 一、杆件在局部坐标系中的刚度矩阵一、杆件在局部坐标系中的刚度矩阵 平面桁架程序计算原理及程序编制 注意注意L= uB-uA，并且由平衡关系得到，并且由平衡关系得到 () ABAB EF FFNuu L 可把可把FA、FB排成列向量排成列向量N，uA， uB，排成列向，排成列向 量量u，系数阵排成矩阵，系数阵排成矩阵RD，即，即 AA BB Fu Nu Fu EFEF LL RD EFEF LL RDuN 平面桁架程序计算原理及程序编制 设杆端分别得到平行于xoy坐标轴的位移uA

20、，vA，uB，vB；则 杆件的伸长量为 cossin () cos() sin BABA Luv uuvv 二、位移的坐标转换 把uA，vA，uB，vB 写成列向量V，系数排列成行向量BT， 上式可以写成如下形式 T LBV 其中 cos sin cos sin B A A B B u v V u v 平面桁架程序计算原理及程序编制 把把L代入上面的公式，得到代入上面的公式，得到 cossincossin A A B B u v EF NRDBV uL v 此时，已用全局坐标系中的位移表达出杆的内力此时，已用全局坐标系中的位移表达出杆的内力N。 平面桁架程序计算原理及程序编制 由于最终得到的平

21、衡方程都是相对于全局坐标系建立的，上面由于最终得到的平衡方程都是相对于全局坐标系建立的，上面 所计算出的内力与杆轴方向一致。假设在杆端作用平行于杆轴所计算出的内力与杆轴方向一致。假设在杆端作用平行于杆轴 的力的力FA、FB，则由平衡方程，则由平衡方程 AB FFN FA和和FB在在x轴、轴、y轴方向的分量分别为：轴方向的分量分别为： coscos AxA FFN sinsin AyA FFN coscos BxB FFN sinsin ByB FFN 平面桁架程序计算原理及程序编制 把把 、 、 、 排列成列向量排列成列向量 ，则有，则有 Ax F Ay F Bx F By F P cos s

22、in cos sin Ax Ay Bx By F F PN F F 由此由此 TT PBNBRDBV 令令 T CBRDB 得到得到 PCV 全局坐标系下的 单元刚度矩阵 平面桁架程序计算原理及程序编制 5-8 结构总刚度矩阵的形成结构总刚度矩阵的形成 假设平衡方程形式为假设平衡方程形式为 R WP 展开上式可以写成 11112211 21122222 1122 nnnn nnnn nnnnnnnnnnnn r wr wr wP r wr wrwP r wrwrwP 式中nn=2*NU；P1，P2，Pnn为加在节点上的外力。 平面桁架程序计算原理及程序编制 5-9 5-9 杆件内力的计算杆件内

23、力的计算 对任意一根杆件，如果杆的两端位移已知。这些位对任意一根杆件，如果杆的两端位移已知。这些位 移在总位移向量中移在总位移向量中i0+1，i0+2，j0+1，j0+2的位置的位置 上，表示如下上，表示如下 0 1 0 2 0 1 0 2 i i j j w w V w w 进而用下面公式计算结构内力进而用下面公式计算结构内力 NRDBV 平面桁架程序计算原理及程序编制 平面桁架计算程序的编制平面桁架计算程序的编制 平面桁架程序计算原理及程序编制 矩阵位移法求解一般步骤：矩阵位移法求解一般步骤： （1 1）结构计算简图。（节点、单元编号，建立一个统一的坐）结构计算简图。（节点、单元编号，建立

24、一个统一的坐 标系等）标系等） （2 2）分析节点位移的力学特性，确定位移未知数。）分析节点位移的力学特性，确定位移未知数。 （3 3）建立每根杆件两端位移和内力的关系。（单根杆件的刚）建立每根杆件两端位移和内力的关系。（单根杆件的刚 度矩阵）度矩阵） （4 4）根据每根杆件上的上述关系建立结构可动节点的平衡方）根据每根杆件上的上述关系建立结构可动节点的平衡方 程。（结构总刚度矩阵）程。（结构总刚度矩阵） （5 5）求解平衡方程，得到节点位移。）求解平衡方程，得到节点位移。 （6 6）根据求得的位移，利用每根杆件位移与内力关系计算各）根据求得的位移，利用每根杆件位移与内力关系计算各 杆的内力。

25、杆的内力。 平面桁架程序计算原理及程序编制 程序流程程序流程 平面桁架计算程序的总功能表示平面桁架计算程序的总功能表示 依据平面桁架及计算模型的有关数据，进依据平面桁架及计算模型的有关数据，进 行结构静力计算，并输出计算结果行结构静力计算，并输出计算结果 平面桁架程序计算原理及程序编制 根据平面桁架的特性进一步展开为：根据平面桁架的特性进一步展开为： 总体数据结构的设计总体数据结构的设计 输入描述平面桁架本体结构和性质的数据输入描述平面桁架本体结构和性质的数据 输出结构图形输出结构图形 计算各杆件的单元刚度矩阵计算各杆件的单元刚度矩阵 形成结构总刚度矩阵形成结构总刚度矩阵 总刚度矩阵三角化总刚

26、度矩阵三角化 输入荷载数据，形成右端项总外力向量输入荷载数据，形成右端项总外力向量 回代求解，求出总位移向量回代求解，求出总位移向量 打印位移打印位移 计算各杆内力计算各杆内力 输出各杆内力输出各杆内力 平面桁架程序计算原理及程序编制 根据平面桁架的特性进一步展开为：根据平面桁架的特性进一步展开为： 总体数据结构的设计总体数据结构的设计 输入描述平面桁架本体结构和性质的数据输入描述平面桁架本体结构和性质的数据 输出结构图形输出结构图形 计算各杆件的单元刚度矩阵计算各杆件的单元刚度矩阵 形成结构总刚度矩阵形成结构总刚度矩阵 总刚度矩阵三角化总刚度矩阵三角化 输入荷载数据，形成右端项总外力向量输入

27、荷载数据，形成右端项总外力向量 回代求解，求出总位移向量回代求解，求出总位移向量 打印位移打印位移 计算各杆内力计算各杆内力 输出各杆内力输出各杆内力 平面桁架程序计算原理及程序编制 总体数据结构的设计总体数据结构的设计 平面桁架程序的总体数据结构设计，包括平面桁架程序的总体数据结构设计，包括 全局量（标识符）的存储安排与说明。全局量（标识符）的存储安排与说明。 平面桁架程序计算原理及程序编制 一个平面桁架结构的数据有三个方面：一个平面桁架结构的数据有三个方面： n 结构本体描述数据（结构本体描述数据（NW, IESG, NU, X, Y, HL, HRNW, IESG, NU, X, Y,

28、HL, HR） n 性质数据（性质数据（EFEF） n 荷载数据（荷载数据（P PX X, P, PY Y） NWNW为节点总数为节点总数 IESGIESG杆件总数杆件总数 NUNU可动节点总数可动节点总数 X, Y X, Y 节点坐标节点坐标 HL, HR HL, HR 每根杆件两端节点编号每根杆件两端节点编号 EF EF 性质数据性质数据 P PX X, P, PY Y 外载荷数据 外载荷数据 平面桁架程序计算原理及程序编制 MODULE PDATAMODULE PDATA REAL, PUBLIC : BAR_HL(100),BAR_HR(100), REAL, PUBLIC : BAR

29、_HL(100),BAR_HR(100), Y X;Y 1 11 1 1 01 0 0 00 0 1 4 100 1 4 100 HI;HL;EF HI;HL;EF 1 2 1001 2 100 2 3 1002 3 100 1 3 1001 3 100 1 0 1 0 PX;PY PX;PY 1 01 0 节点总数：节点总数： NW=4NW=4 可动节点数：可动节点数： NU=2NU=2 位移未知数总数：位移未知数总数：NDISP=2NDISP=2* *NU=4NU=4 平面桁架结构平面桁架结构 1 34 2 平面桁架程序计算原理及程序编制 输入描述平面桁架本体结构和性质的数据输入描述平面桁

30、架本体结构和性质的数据 结构本体描述数据（结构本体描述数据（NW, IESG, NU, X, Y, HL, HR） 性质数据（性质数据（EF） 平面桁架程序计算原理及程序编制 READ(1,READ(1,* *) NW) NW READ(1,READ(1,* *) IESG) IESG READ(1,READ(1,* *) NU) NU NDISP=NU+NUNDISP=NU+NU WRITE(2,WRITE(2,* *) NW=,NW, IESG=,IESG, NU=,NU) NW=,NW, IESG=,IESG, NU=,NU WRITE(2,WRITE(2,* *) NODE X Y)

31、NODE X Y DO 10 I=1,NWDO 10 I=1,NW READ(1, READ(1,* *) NODP_X(I),NODP_Y(I) NODP_X(I),NODP_Y(I) WRITE(2, WRITE(2,* *) I,NODP_X(I),NODP_Y(I) I,NODP_X(I),NODP_Y(I) CONTINUECONTINUE WRITE(2,WRITE(2,* *) ELEMENT HL HR EF) ELEMENT HL HR EF DO 20 NT=1,IESGDO 20 NT=1,IESG READ(1, READ(1,* *) BAR_HL(NT),BAR_H

32、R(NT),BAR_EF(NT) BAR_HL(NT),BAR_HR(NT),BAR_EF(NT) WRITE(2, WRITE(2,* *) BAR_HL(NT),BAR_HR(NT),BAR_EF(NT) BAR_HL(NT),BAR_HR(NT),BAR_EF(NT) 20 CONTINUE ! LIST DATA, NOW TO ORGANIZE THE GLOBAL MATRIX20 CONTINUE ! LIST DATA, NOW TO ORGANIZE THE GLOBAL MATRIX 输入本体结构和性质的数据输入本体结构和性质的数据数据数据 平面桁架程序计算原理及程序编制

33、4 4 NW NW 4 4 IESG IESG 2 2 NU NU 0 1 0 1 X;Y X;Y 1 11 1 1 01 0 0 00 0 1 4 100 1 4 100 HI;HL;EF HI;HL;EF 1 2 1001 2 100 2 3 1002 3 100 1 3 1001 3 100 1 0 1 0 PX;PY PX;PY 1 01 0 平面桁架程序计算原理及程序编制 根据平面桁架的特性进一步展开为：根据平面桁架的特性进一步展开为： 总体数据结构的设计总体数据结构的设计 输入描述平面桁架本体结构和性质的数据输入描述平面桁架本体结构和性质的数据 输出结构图形输出结构图形 计算各杆件

34、的单元刚度矩阵计算各杆件的单元刚度矩阵 形成结构总刚度矩阵形成结构总刚度矩阵 总刚度矩阵三角化总刚度矩阵三角化 输入荷载数据，形成右端项总外力向量输入荷载数据，形成右端项总外力向量 回代求解，求出总位移向量回代求解，求出总位移向量 打印位移打印位移 计算各杆内力计算各杆内力 输出各杆内力输出各杆内力 平面桁架程序计算原理及程序编制 计算各杆件的单元刚度矩阵计算各杆件的单元刚度矩阵 EF RD L T CBRDB cos sin cos sin B A A B B u v V u v 相应位移向量相应位移向量 平面桁架程序计算原理及程序编制 SUBROUTINE STIF(NT,L,C,B,RD

35、) USE PDATA INTEGER NT,I,J REAL L,RD,C(4,4),B(4),CA,SA,DX,DY I=BAR_HL(NT) J=BAR_HR(NT) !LEFT AND RIGHT NODE NUMBER DX=NODP_X(J)-NODP_X(I) DY=NODP_Y(J)-NODP_Y(I) L=SQRT(DX*DX+DY*DY) CA=DX/L; SA=DY/L !(* DIRECTIN COSINE *) B(1)=-CA B(2)=-SA B(3)=CA B(4)=SA ! (* THE DISPLACEMENT-DEFORMATION TRANSFORMAT

36、ION MATRIX *) RD=BAR_EF(NT)/L !(* BAR STIFFNESS*) DO 10 I=1,4 DO 10 J=1,4 C(I,J)=RD*B(I)*B(J) CONTINUE ! (* C IS THE STIFFNESS MATRIX IN GLOBAL COORDINATE *) RETURN END ! (*STIF*) 计算各杆件的单元刚度矩阵计算各杆件的单元刚度矩阵 平面桁架程序计算原理及程序编制 根据平面桁架的特性进一步展开为：根据平面桁架的特性进一步展开为： 总体数据结构的设计总体数据结构的设计 输入描述平面桁架本体结构和性质的数据输入描述平面桁架本

37、体结构和性质的数据 输出结构图形输出结构图形 计算各杆件的单元刚度矩阵计算各杆件的单元刚度矩阵 形成结构总刚度矩阵形成结构总刚度矩阵 总刚度矩阵三角化总刚度矩阵三角化 输入荷载数据，形成右端项总外力向量输入荷载数据，形成右端项总外力向量 回代求解，求出总位移向量回代求解，求出总位移向量 打印位移打印位移 计算各杆内力计算各杆内力 输出各杆内力输出各杆内力 平面桁架程序计算原理及程序编制 形成结构总刚度矩阵形成结构总刚度矩阵 对总刚度矩阵对总刚度矩阵R清零清零 FOR NT:=1IESG DO 注：对所有杆件循环注：对所有杆件循环 BEGIN 调出调出STIF(NT,T,C,B,RD)；计算杆件

38、在全局坐标系下单元刚度；计算杆件在全局坐标系下单元刚度 矩阵矩阵C 调出调出I0J0(NT,II)；计算对号入坐信息数组；计算对号入坐信息数组II 调出调出ASSEMB2(IB,C,R,NDISP,BV,V)； 注：结构总刚度矩阵注：结构总刚度矩阵R 已形成已形成 END; 平面桁架程序计算原理及程序编制 ! !确定这根杆件两端点（确定这根杆件两端点（HLHL端和端和 HRHR端）位移在总位移向量中的编号端）位移在总位移向量中的编号 SUBROUTINE I0J0(NT,II) USE PDATA INTEGER NT,II(4),I,I0,J0 I=BAR_HL(NT)-1 I0=I+I !

39、 NOTE THE NUMBERING FORMULA DO 10 I=1,2 10 II(I)=I0+I I=BAR_HR(NT)-1 J0=I+I DO 20 I=3,4 II(I)=J0+I-2 RETURN END ! ( *I0J0 *) uiui在总位移向量中是第在总位移向量中是第2 2* *（i-1i-1）+1+1号号 ujuj在总位移向量中是第在总位移向量中是第2 2* *（j-1j-1）+1+1号号 平面桁架程序计算原理及程序编制 SUBROUTINE STIF(NT,L,C,B,RD) USE PDATA INTEGER NT,I,J REAL L,RD,C(4,4),B(

40、4),CA,SA,DX,DY I=BAR_HL(NT) J=BAR_HR(NT) !LEFT AND RIGHT NODE NUMBER DX=NODP_X(J)-NODP_X(I) DY=NODP_Y(J)-NODP_Y(I) L=SQRT(DX*DX+DY*DY) CA=DX/L; SA=DY/L !(* DIRECTIN COSINE *) B(1)=-CA B(2)=-SA B(3)=CA B(4)=SA ! (* THE DISPLACEMENT-DEFORMATION TRANSFORMATION MATRIX *) RD=BAR_EF(NT)/L !(* BAR STIFFNES

41、S*) DO 10 I=1,4 DO 10 J=1,4 C(I,J)=RD*B(I)*B(J) CONTINUE ! (* C IS THE STIFFNESS MATRIX IN GLOBAL COORDINATE *) RETURN END ! (*STIF*) 平面桁架程序计算原理及程序编制 1358101118222533 0112106327 对角元地址数组： 半带宽数组： 3.2 -1.04.3 0.82.0 4.00.18.8 4.14.0 2.1 -0.81.20.30002.8 1.5003.1 0.703.7 -0.200-1.20-2.14.17.2 刚度矩阵（二维存储）

42、： 3.24.30.82.04.00.18.84.14.0 2. 1 1.20.300002.81.5003.10.703.70004.17.2 刚度矩阵（一维存储）： 平面桁架程序计算原理及程序编制 存储每一行的对角元素在一维存储数组中存储每一行的对角元素在一维存储数组中R中的位中的位 置，存放在数组置，存放在数组 VBV:BV+NDISP 对角元素在一维数组对角元素在一维数组R中的地址为中的地址为VBV+I，而这，而这 个元素的值为个元素的值为 RVBV+I 1358101119232634 对角元地址数组：对角元地址数组： 3.24.30.82.04.00.18.84.14.0 2. 1

43、 1.20.300002.81.5003.10.703.70004.17.2 一维刚度数组：一维刚度数组： 平面桁架程序计算原理及程序编制 根据根据V可以求出：可以求出： （1）第）第i行主元在行主元在R中的位置为中的位置为VBV+i，数值为，数值为RVBV+i。 如第如第3行，行，BV=0，i=3，V3=6，其值，其值R6=r33。 （2）第）第i行半带宽（不计对角元）为：行半带宽（不计对角元）为： BDWi= VBV+i- VBV+i-1-1 （3）在方阵中处于）在方阵中处于i行行j列的元素列的元素rij在在R中的地址为：中的地址为： ix= VBV+i-i+j 当当i=j 但是要求但是要

44、求ix VBV+i-1 （i=j），否则该元素为零，在），否则该元素为零，在R中并没有中并没有 存放。存放。 （4）第）第i行第一个非零元素所在列号为行第一个非零元素所在列号为 i1= VBV+i-1- VBV+i+i+1 式中当式中当i=1时时VBV定义为零。定义为零。 （5）第）第i行对角元地址可用下面的递推公式计算：行对角元地址可用下面的递推公式计算： Vi=Vi-1+BDWi+1 式中式中BDWi为第为第i行的半带宽。行的半带宽。 平面桁架程序计算原理及程序编制 半带宽选大的自然语言表示如下：半带宽选大的自然语言表示如下： 未知数总数未知数总数NDISPNDISP已知已知 安排数组安排

45、数组V V的存区的存区 半带宽存区清零，准备选大半带宽存区清零，准备选大 对杆件循环对杆件循环 调出调出I0J0I0J0（NT, IINT, II） 为计算为计算NTNT号杆对半带宽作准备号杆对半带宽作准备 选出半带宽选出半带宽 平面桁架程序计算原理及程序编制 根据自然语言表示，半带宽选大的程序段是根据自然语言表示，半带宽选大的程序段是: : BV:=1;BV:=1; FOR I:=0 TO NDISP DOFOR I:=0 TO NDISP DO VBV+I:=0VBV+I:=0 注：半带宽选大之前，对注：半带宽选大之前，对V V数组清零数组清零 FOR NT:=1 TO IESG DOFO

46、R NT:=1 TO IESG DO BEGINBEGIN I0J0(NT,II)I0J0(NT,II) 杆位移对号杆位移对号 FOR I:=1 TO IB DO ! IBFOR I:=1 TO IB DO ! IB为数组为数组C C的下标界偶，即单元位移向量的元素数的下标界偶，即单元位移向量的元素数 BEGIN FOR J:=1 TO I DOBEGIN FOR J:=1 TO I DO 对杆件的行循环和列循环对杆件的行循环和列循环 IA:=III;IA:=III; JA:=IIJ;JA:=IIJ; 转换到总刚度矩阵的行、列号转换到总刚度矩阵的行、列号IAIA与与JAJA IF IAJA T

47、HENIF IA0) AND (IA=NDISP) THEN 防止出界防止出界 BEGIN IF VBV+IAIA-JA THEN VBV+IA:=IA-JA; ! 选大选大 END; END; END; 半带宽已选出半带宽已选出 平面桁架程序计算原理及程序编制 在此基础上就可以编制在此基础上就可以编制 DIAGADR(NDISP,V) 形成总刚度矩阵的对角元地址过程。形成总刚度矩阵的对角元地址过程。 功能：在给定总刚度矩阵功能：在给定总刚度矩阵R存放的基址存放的基址IV，半带宽数组，半带宽数组 VBV:BV+NDISP的条件下，求出总刚度矩阵各行的对角元的条件下，求出总刚度矩阵各行的对角元

48、地址，并仍就存放在地址，并仍就存放在VBV+1 VBV+NDISP中。中。 平面桁架程序计算原理及程序编制 SUBROUTINE ASSEMB2(IB,II,C,R,NDISP,V,PHASE) USE PDATA INTEGER IB,NDISP,PHASE,I,J,K,IA,JA,II(4),V(120) REAL C(4,4),R(1200) DO 100 I=1,IB DO 50 J=1,I IA=II(I) JA=II(J) ! TAKE OUT DISP NUMBER IF (IA .LT. JA) THEN K=IA IA=JA JA=K ENDIF IF (JA.GT.0).A

49、ND.(IA.LE.NDISP) THEN IF (PHASE.EQ.2) THEN K=V(IA)-IA+JA ! ( *THE ADESS IN R * ) R(K)=R(K)+C(I,J) ! ( *ACCUMULATING * ) ELSE IF (V(IA).LT.(IA-JA) V(IA)=IA-JA END IF END IF 50 CONTINUE 100 CONTINUE RETURN END ! ASSEMB2 平面桁架程序计算原理及程序编制 SUBROUTINE DIAGADR(NDISP,V) USE PDATA INTEGER V(120),IA ! FROM SEM

50、I-BANDWIDTH TO DIAGONAL ADRESS OF GLOBAL MATRIX ! IV - - BASE ADRESS OF GLOBAL MATRIX ! VBV.BV+NDISP SEMI-BAND IN BUT RETURN WITH DIAGONAL ADRESS ! NDISP : GLOBAL VARIABLE DO 10 IA=1,NDISP 10 V(IA)=V(IA)+V(IA-1)+1 RETURN END ! DIAGADR 平面桁架程序计算原理及程序编制 总刚度矩阵三角化总刚度矩阵三角化 修改的平方根法运算工作量小且可以处理满阵修改的平方根法运算工作量

51、小且可以处理满阵 存储的系数矩阵，可以使用较小的存储空间。存储的系数矩阵，可以使用较小的存储空间。 每增加一个右端项时只须增加一次回代求解，每增加一个右端项时只须增加一次回代求解， 而无须重新分解系数矩阵。具体实施在形成结而无须重新分解系数矩阵。具体实施在形成结 构总刚度矩阵之后。构总刚度矩阵之后。 ! NDISPTHE ORDER OF MATRIX TO BE LDLTED ! VBV.BV+NDISPTHE ARRAY OF DIAGONAL ADREDD, THE MATRIX TO BE ! LDLTED IS IN RVBV+1.VBV+NDISP; ! T1.BDWTHE WORKING ARRAY; BDW IS THE MAX-BAND OF R SUBROUTINE LDLT1(NDISP,V,R,T) 平面桁架程序计算原理及程序编制 根据平面桁架的特性进一步展开为：根据平面桁