第二部分 平面杆系结构的线弹性有限元法.doc

3 平面杆系结构的线弹性有限元法3.1概论在有限元法中，可以采用位移法，也可以采用力法或混合法。其中提出最早并且应用最广的是位移法。对于平面杆系结构来说，位移法实际上就是结构力学中的矩阵位移法（也称刚度法），在计算时以结点位移作为基本未知量。杆系结构的矩阵分析实际上就是有限元法。其基本思路是：先把结构离散成有限个数目的单元，然后再考虑某些条件，将这些离散的单元重新组合在一起进行分析计算。这样使一个复杂的计算问题转化为简单的单元分析和集合问题。根据这个思路，杆系结构的有限元法可分为两大步骤：（1）单元分析。研究单元的受力与变形之间的关系；（2）整体分析。研究如何将这些离散的单元重新组合得到与实际问题相符合的（如边界条件、外界荷载等等）的计算模型整体刚度方程。在有限元中，一般采用矩阵形式进行分析求解，因为矩阵运算不仅使公式非常紧骤，而且形式统一，易于编程，适合在电子计算机上进行自动求解。因此，在有限元法的一般格式中，应尽量采用矩阵形式进行运算。3.2 局部坐标系下的单元刚度矩阵1单元的划分。在杆系结构的有限元法中，一般将由相同材料、具有相同横截面的一根杆件（即等截面直杆）当成一个单元，整个结构就是由有限个杆件单元组成的集合体。杆件单元具有2个结点，即首结点和末结点，但一般是先确定结点的位置，结点一旦确定，则结点之间的单元也就确定了。在进行杆系结构的单元划分时，应注意如下事项：结点位置的确定。结点一般选在杆件的如下位置：杆件的转折点、杆件汇交点、支承点、截面或材料的突变点，这些点都是结构的构造点，有时为了使结构只承受结点荷载，在集中荷载的作用处也设置一个结点。结点的编号。为了使集合以后的总刚的带宽最小，一般应遵循尽量使相关结点（有单元相连的结点）编号差值的最大值最小的原则进行。2单元刚度矩阵图1 梁单元考虑一等截面的平面梁单元，单元首末结点分别为，单元长为，单元抗弯刚度为，为材料的弹性模量，是截面的抗弯惯矩，取轴为沿梁单元中心轴，轴与轴成90o，如图1所示。位移模式和形函数如果不考虑杆件的轴向变形与横向弯曲变形的互相影响，且设轴向的位移（即单元轴向位移）取为的线性函数，而对于轴向的位移（即单元横向位移，亦即梁单元挠度）取为的三次多项式表示。于是有： (1)不考虑剪切变形的影响，即假定，可推导出位移模式又可表示为（2）其中，其中；式中、反映了单元的位移形态，称为形态函数，简称形函数，称为形函数矩阵。用结点位移表示的应变和应力如单元产生拉压变形和弯曲变形，则其纵向纤维的线应变可分成两个部分：（1）称为拉压应变，也称为轴向应变，整个梁截面都相同；（2）弯曲应变，沿梁高因存在曲率（称为广义弯曲应变）而不同。假设不考虑剪切变形，由梁单元的几何方程可知（3）在杆件截面高度处，有曲率引起的应变为 （4）由（3）和（4）式可写成（5）式中表示对求一阶导，表示对求二阶导，其余说明类似。上式可写成（6）由虎克定律，就可以得到利用结点位移表示的应力表达式 （7）式中称为弹性矩阵，且由式（6）可得（8）式中，称为应力矩阵。由虚位移原理导出梁单元的刚度方程应用虚功原理于梁单元上，可得到梁单元保持平衡的刚度方程为（9）令 （10） （11）式中称为单元刚度矩阵，为由于分布荷载而移置的等效结点荷载，为直接作用在结点上的荷载。于是 （12）令（13）则（14）上式即为局部坐标系中的单元刚度方程。将的具体的表达式代入（12）式，经积分和矩阵运算可得到平面梁单元的单元刚度矩阵为,具体表达式如下:（15）式(13)中元素为平面杆系单元刚度矩阵元素，即：显然单元刚度矩阵为对称矩阵。3.3 整体坐标系下的单元刚度矩阵在前面的分析之中，单元刚度矩阵是在单元的局部坐标系中形成的，由于各个单元的局部坐标系不同，因此必须将每个单元的单刚转换到同一个公共的坐标系下，这个公共坐标系就是整体坐标系。为了区别起见，在局部坐标系下的杆端分量符号顶上加“”。下面首先介绍转换矩阵的概念，然后据之建立整体坐标系下的单元刚度矩阵。一转换矩阵如图32所示，任一单元的首端结点力在两种坐标系中的分量。其中图(a)表示在局部坐标系中的三个分量、和，图(b)表示在整体坐标系中的三个分量、和。为了导出、与、之间的关系式，在图中将两个力、分别投影在和轴上，可得出下式中的前两式：其中第三式表明，在两个坐标系中的力偶仅是仍彼此相等。表示由轴转到轴的角度，以逆时针方向为正。同理对单元的另一端力也可得出类似的关系：把以上两个方程组组合成一个矩阵方程得或简写成式中是局部坐标系中的单元杆端力列阵，是整体坐标系中的单元杆端力列阵，T称为单元的坐标转换矩阵。可以证明，T矩阵是正交矩阵，其逆矩阵等于它的转置矩阵，即对于单元杆端结点位移，也可以同样进行转换，即有二整体坐标系中的单元刚度矩阵1推导过程。现在研究整体坐标系中的单元等效杆端力与单元结点位移之间的关系式。因为所以有令则上式即为整体坐标系下的单元刚度方程，简称单元整体刚度方程。因此相应地称为单元整体刚度矩阵。2单元刚度矩阵的特点。（1）由，可知单元刚度矩阵为对称矩阵。（2）单元刚度的分块由于在以后的整体分析之中，是对结构的每个结点建立平衡方程，为了以后讨论方便，可把上式按单元的首未结点编号i、j进行分块，写成如下形式：式中：分别为首端i和末端j的杆端力和杆端位移列阵。为单元刚度矩阵的四个子块，每个子块都是阶的方阵，并且有，两个子块的元素相同，子块为对称子块。同时可以得到如下式子：3.4整体刚度矩阵（即总刚）一回顾在结构力学中的位移法一章中，我们建立位移法的典型方程是通过附加力臂的力矩平衡或杆件在某一方向的受力平衡来得到的。平面刚架的整体刚度方程也是通过结点的受力平衡得到的。二推导假设先不考虑结构的支承条件和约束条件，结构只承受结点荷载的情况。设有刚结点i的相关结点为j、k、l，则相应的单元为命名为ij、ik、il，如图所示：设结点i承受的结点荷载为：在平衡状态下，结点i可以建立三个平衡方程，即：上式可以写成矩阵向量形式，有将与I相关的单元的杆端力向量代入上式可以得到：很显然，有故有亦即上式包含三个方程，。N为结构总结点数。由此我们可以知道，如果一个结构由N个结点组成，则可以建立3N个方程。由这3N个方程组成的方程组的系数矩阵称为整体刚度矩阵，简称总刚，用来表示。该矩阵按上式可以分块为N行N列个子块组成。上式中各项系数子块就是组成刚度矩阵相应的第i行j、k、l列的子块元素（这里的行和列是以子块作为一个元素来划分的，实际上，每个子块还有3行3列）。是i行i列主对角线上的子块，实际上是i行j、k、l列上的子块。如果i结点再没有其它相关结点。则该行其它列上的子块为零子块。由上述分析，我们可以初步得到总体刚度矩阵形成的规律。为了更形象地说明，下面我们采用如下示意图来说明之。单元ij是由首结点i、末结点j组成（首结点编号不一定少于未结点编号，且先假定ij），图中i03(i-1)、j03(j-1)，单刚按照“对号入座”的原则，由箭头示意置于总刚中的相应位置。同理，单元ik也可以按如图所示的“对号入座”累加到相应的总刚中的位置。为了讨论方便，将主对角线上的子块称为主子块，其余的子块称为副子块，同交于一个结点的各杆单元称为该结点的相关单元。三结论由上面分析和说明，我们可以得到如下的结论：（1）、总刚中的主子块是由结点i的各相关单元的主子块累加而成，亦即。（2）、总刚中的副子块，当i、m为相关结点时，等于im单元的相应的副子块，即，如果i、m不为相关结点，即没有单元相连，则为零子块。（3）、总刚的形成可以按如下方式形成：将所有单元的单刚的相应子块按上述“对号入座”的原则累加到相应的总刚位置上，全部累加完成，总刚也就形成。在编写程序时，用循环形式极容易完成。四总刚的特点(1)、是对称矩阵，因为单刚的子块，故有总刚中的相应副子块与关于对角线对称，为对称子块；主子块，由前面分析可知，因为为对称子块，故亦为对称子块，所以整个刚度矩阵为对称矩阵。(2)、是奇异矩阵，因为建立整体刚度方程的过程，我们没有考虑结构与基础的连接情况。也就说，这样分析的结构是没有支承点的，具有无穷多个位移解。因而肯定是奇异矩阵，由于此时总刚中没有引入边界条件，故有时也称之为原始刚度矩阵。(3)、是带状矩阵。由前面总刚形成的规律可以看出，第i行（以单刚中的子块为一个元素）中非零子块的最大列号或最小列号m（当mi时，为最大；当mi时为最小）是i结点的一个相关结点的编号，很显然，当i、m编号相差最小时，带宽最小。因此，为了节省存储空间，应尽量使结构中相关结点编号的差值最小，这是杆系结构进行有限元计算时结点编号的一个总原则。例4-1：见结构力学（李廉锟编）P249页。试求图示刚架的原始总刚。各杆材料和截面均相同，解：（1）、将各单元、结点编号，并选取整体坐标系和各单元的局部坐标系如图所示。（2）、各单元的整体坐标系中的单刚按前面公式计算，先计算所需相关数据如下：对于单元，可计算得对于单元和，可计算得（3）、将以上各单刚子块对号入座即得总刚：3.5 支承条件的引入和非结点荷载的处理一 结构支承条件的引入在上述分析过程中，没有考虑支承条件和约束情况。实际上，支承条件也有在总刚形成之前就考虑进去，这是先处理法的思路。采用这种方法时，在准备工作阶段，不仅要对结点、单元进行编号，而且还要对结点位移未知量统一编号。不发生位移的结点自由度方向的位移编号为0，不作为结点未知量来考虑。采用这种方法编写程序稍微复杂一些。为了与前面讲述的内容一致，我们采用后处理法。所谓后处理法，就是在原始刚度矩阵形成之后，根据结构的支承条件，对原始总刚和荷载列阵进行处理，从而得到引入了支承条件的刚度方程。下面根据后处理法讲述支承条件的处理方法。1支承输入信息的表达与支座类型：我们知道，在结构力学中，结构与地基的连接装置可以简化成各种支座，由于地基的不均匀沉陷，引起各种支座强迫位移。因此，必须将结构的全部约束信息通过某种形式传递给计算机。一般在结点编号时，在杆件与支座相连处设置成一个结点，此结点称为支座结点。首先，必须明白结构有多少个结点受到了约束，受约束的结点编号。其次，每个受约束的结点在哪个自由度（即方向的线位移及转角，有时简称结点的1、2、3自由度）受到约束，是否有强迫位移，若有，其值为多少。如果支座结点的某自由度方向的位移为零，则该方向的约束信息填0.0，如果没有约束，则置一大数9999.0（因为结构符合线弹性、小变形假设，一般结点位移小于此值），如果有支座强迫位移，则填具体的位移值。如铰支座结点，其，则结点约束信息填0.0、0.0、9999.0。下面为右图所示刚架支承约束信息的一种输入方式：2 （结构受约束的结点个数）10.00.00.0（第一个约束结点的支承信息）40.00.00.0（第二个约束结点的支承信息）：常见的各种支座形式和约束信息见下表。支座名称简 图约束信息UV固定支座0.00.00.0铰支座0.00.09999滚轴支座199990.0999920.099999999滑动支座199990.00.020.099990.02处理办法边界条件的处理方法通常有如下几种：（1）缩减总刚法。对于支承结点某个自由度方向位移为零时，划去总刚中与该自由度相对应的行和列，同时划去荷载列阵中相应的行。这样就对总刚和荷载列阵进行了阶次缩减，降低了方程的阶次，从而提高了计算速度。但紧缩总刚和荷载列阵的程序较为复杂。（2）充0置1法。对于结构某个支承自由度方向的位移为已知（通用于零位移和发生支座强迫位移两种情况）。为了使方程求解出来的该自由度i方向的位移为，即，需要将荷载列阵中用取代，整个刚度方程如下式所示：（3）、乘大数法。如果结构在某个支承自由度方向的位移为已知（通用于零位移和发生支座强迫位移两种情况）时，同样采用乘大数的方法。具体做法是将第i个方程进行如下处理：上式中除了包含M的两项外，其它各项相对于它们都比较小，可以忽略不计。因此，上式即归结为给定的支承条件。（4）、作为非结点荷载处理。把支座结点位移转换成与该结点相连的各单元的杆端位移，将单元看成两端固定梁，在局部坐标系中求出单元在给定杆端位移下的固端力，求出单元的等效结点荷载，然后计算。二 非结点荷载的处理所谓非结点荷载就是指作用在杆单元中间的集中力或分布力。1思路在实际问题中，不可避免地遇到非结点荷载。在连续梁一章的学习中我们已经知道，对于这种情况，可以利用叠加原理，分两步进行处理：（1）与位移法类似，首先在结构上加附加刚臂和附加链杆，阻止所有结点的线位移和角位移，此时，每个杆件单元都成了两端固定的固支梁，各单元两端承受固端力，附加约束上也施加了反力和反力矩。反力和反力矩的大小可以根据结点平衡求得，其大小就等于汇交于该结点的各固端力的代数和。（2）取消附加约束。取消附加约束的方法是在各附加约束上加上与上述反力和反力矩大小相等、方向相反的力和力矩，此力和力矩就称为原非结点荷载的等效结点荷载。（3）前面两种状态的叠加，结构的内力和变形与原来的实际情况是相同的。因此，内力和变形也相应地由（1）、（2）两种情况下的结果叠加而成。第一种情况下，各杆件单元的杆端内力就是固端力，一般可以在结构力学教材书上查阅而得，无需计算。第二种情况下各杆单元的内力就是结构在结点荷载（包括等效结点荷载和直接作用在结点上的荷载）单独作用下的内力，需要通过有限元计算得到。对应于（1）、（2）、（3）三种受力状态，可用如下结构图a、b、c来说明。图a所示刚架的受力和变形可由图b和图c两种状态的叠加而成。图b所示的状态其内力图可直接查表可得，其中所有单元的杆端力其实就是固端力。仅仅需要计算图c所示状态下的结构内力，该状态下结构所受的结点荷载与图b中各附加约束的反力和反力矩大小相等、方向相反。此时内力需通过计算得到。这些结点荷载就是图a中杆件单元非结点荷载的等效结点荷载。2等效结点荷载列阵的推导根据前面的分析和图解说明，设某单元e的首尾结点编号为i、j，杆中承受非结点荷载，在局部坐标系中其固端力为F表示固端的意思（“Fixed end”的首个字母），这些固端力很容易由公式或表格查得。将局部坐标系下的固端力转换到整体坐标系下，则有：任一结点i的等效结点荷载等于汇交于这一结点的各杆件单元的固端力的代数和的相反数。如果用表示该结点的等效结点荷载（英文单词“Equivalent”首个字母），即：如果结点i还承受了直接作用在结点上的荷载，其表达式为：则结点i承受的总荷载为：称为结点的综合结点荷载，所有结点的综合结点荷载按顺序形成荷载列阵，则整个结构的综合结点荷载列阵可写成为：3非结点荷载的类型略4结构内力的计算由前面的分析可知，单元e杆端内力（全部转换到局部坐标系下）由两部分组成：一部分是所有结点位移全部约束时的固端力。另一部分是综合荷载作用下单元杆端位移引起的杆端力。公式如下：（1）第一部分，固端力为：（2）第二部分，杆端位移引起的杆端力为：（3）杆件单元e的最终杆端内力为： 3.6矩阵位移法的计算步骤和示例通过前面的分析和讨论，我们可以归纳出利用矩阵位移法进行结构分析的步骤：（1）准备数据。对结构的结点和单元编号（结点编号时遵循相关结点编号差值的最大值尽量少）；选定整体坐标系和局部坐标系；确定各结点的x、y坐标值；计算各种类型的材料参数（弹性模量、容重等等）和几何截面参数（截面积、抗弯惯性矩和抗弯截面模量等等）；确定各单元的首尾结点编号以及所属的材料类型号和几何参数类型号；确定各个非结点荷载所作用的单元编号、所属的荷载类型号、荷载大小和位置（以方便等效结点荷载的计算），确定各个直接结点荷载所作用的结点编号、大小和方向。（2）计算各杆件单元的整体坐标系下的单元刚度矩阵。（3）按单元编号循环累加单刚形成原始刚度矩阵。（4）计算固端力、等效结点荷载和综合结点荷载。（5）引入支承条件，修改结构的原始刚度方程。（6）求解刚度方程，求解结点位移。（7）计算各单元在局部坐标系下的杆端力（注意须叠加上各杆件单元的固端力）。例4-2 试求图示刚架的内力，各杆件的材料和截面相同，具体数据见例4-1。解：（1）将结点、单元进行编号，确定坐标系，如图示。（2）求出各单元在整体坐标系下的单元刚度矩阵，见例4-1。（3）将各单刚子块对号入座，形成结构的原始刚度矩阵，见例4-1。（4）计算非结点荷载作用下的各单元固端力、等效结点荷载和综合结点荷载。根据结构力学教材书上的表格可查得，各单元在其局部坐标系中的固端力为：将单元的，单元、的代入计算，可得各单元在整体坐标系下的固端力列向量：由上式可以求出结点2、3上的等效结点荷载为：综合结点荷载为：于是结构的结点外力列向量为：结构的原始刚度方程为：由于现在采用的是手解方程，所以采用缩减矩阵法引入支承条件，支承条件为：在原始刚度矩阵中删除与已知零位移相对应的行和列，同时在结点位移列向量和结点荷载列向量中删除相应的行和列，原始刚度方程修改如下：（6）、解方程，求得结点位移为：（7）、计算各单元杆端力：单元:单元:单元:最后刚架的变形图和弯矩图如下：第三节补充说明3.7平面刚架程序框图一总框图与程序标识符：1总框图开 始结 束输入基本数据0输出计算结果6引入边界条件3计算结构的内力5求解刚度方程4形成荷载列向量2形成总刚矩阵1结 束平面刚架的整个计算过程见上面总框图。以上为初级程序框图。每个初级程序框图还可以进一步分成更细的几个二级程序框图。按照上述总框图编制的平面刚架静力计算程序的适用范围是：（1）、结构形式由等截面直杆组成的具有任意几何形状的平面杆系结构：刚架、组合结构、桁架、排架和连续梁。平面刚架单元之间的连接结点可以是刚结点、铰结点和刚铰混合结点。（2）支座形式结构的支座可以是固定支座、铰支座、滚轴支座和滑动支座。（3）、荷载类型作用在结构上的荷载包括结点荷载和非结点荷载，各种非结点荷载类型见表。（4）、材料性质结构的各个杆件可以用不同的弹性材料组成。在平面刚架的矩阵分析中，考虑了杆件的弯曲变形和轴向变形，而忽略了剪切变形的影响。2程序标识符现将子框图和源程序中主要标识符的意义说明如下：（1）、整型变量NJ结点总数；NM单元总数；N结点位移未知量总数，即整个结构的总自由度数，亦即总刚矩阵的阶数；MB最大的半带宽；NME材料种类数；NGE截面几何特性种类数；NS受约束的结点位移数；NSJ支承结点数；NJP结点荷载个数；JN荷载作用的结点号；JD荷载作用的方向号（1水平力；2垂直力；3集中力偶）；NLM非结点荷载个数；ILT非结点荷载类型号；K单元序号；（2）、整型数组ISE（NM，2）各单元首端和未端结点号；IMG（NM，2）各单元的材料特性号和截面几何特性号；IS（NS）各约束自由度的总体编号；（3）、实型数组WE（NME，2）各种材料特性的容重和弹性模量；AI（NGE，2）各种截面几何特性的截面积和抗弯惯性矩；XY（NJ，2）各结点坐标；SD（NS）各约束自由度的约束位移值（已知的支座位移）；P（N）结构总体荷载列阵，总刚求解后为结点位移增量；TF（NM，6）各单元杆端力增量；BP（NS）各约束点的约束反力增量；DT（N）结点位移累计量；TFT（NM，6）各单元杆端力累计量；二子框图下面讨论主要的子框图。1总刚的形成单刚累加形成总刚矩阵形成局部座标系下的单刚，并按公式得到整体座标系下单元刚度矩阵计算座标转换矩阵计算单刚中的K1K4K1，NM计算单元的几何特性3、子程序功能：计算单元长度、局部坐标系与整体坐标系夹角正弦值和余弦值Subroutine ch1(k,I0,J0,SI,CO,RL)Common /m1/ise(500,2),xy(500,2)调用单元的首尾结点编号计算单元杆件的余弦和正弦计算单元的长度I=ise(k,1)J=ise(k,2)I0=3*(I-1)J0=3*(J-1)Co=xy(J,1)-xy(I,1)Si=xy(J,2)-xy(I,2)Rl=sqrt(co*co+si*si)If(Rl.lt.1.0e-10) write(*,10) k,co,si10format (1x,ie=,i3,5x,dx=,e12.4,dy=,e12.4)Co=co/RlSi=si/Rlend子程序：功能：计算局部坐标系下单刚的4个独立元素。subroutine dke(k,nme,nge,rl,ai,we)common /m2/img(500,2),tf(500,6),ek(500,4)调用单元的几何类型号与材料类型号 dimension ai(nge,3),we(nme,2)m=img(k,2)if(m.le.0) goto 100a=ai(m,1)e=we(img(k,1),2)r=ai(m,2)根据公式计算独立的四个元素ek(k,1)=e*a/rlek(k,2)=12.0*e*r/(rl*rl*rl)ek(k,3)=6.0*e*r/(rl*rl)ek(k,4)=4.0*e*r/rl100continue end4、6、子程序功能：由局部坐标系下的单刚形成整体坐标系下的单刚。subroutine ke(k,t,ae)将局部坐标系单刚所有的的元素置将单刚的上三角中不为0 的元素赋值对单刚的下三角元素赋值将单刚转换成整体坐标系下的单刚common /m2/img(500,2),tf(500,6),ek(500,4)dimension ae(6,6),t(6,6),t2(6,6)call clear(6,6,ae)ae(1,1)=-ek(k,1)ae(1,4)=-ek(k,1)ae(2,2)=ek(k,2)ae(2,3)=-ek(k,3)ae(2,5)=-ek(k,2)ae(2,6)=ek(k,3)ae(3,3)=ek(k,4)ae(3,5)=-ek(k,3)ae(3,6)=ek(k,4)/2.0ae(4,4)=ek(k,1)ae(5,5)=ek(k,2)ae(5,6)=-ek(k,3)ae(6,6)=ek(k,4)do 30 I=1,5do 30 j=I+1,630 ae(j,I)=ae(I,j) call matmul(6,6,6,ae,t,t1) do 20 I=1,6do 20 j=1,620 t2(I,j)=t(j,I) call matmul(6,6,6,t2,t1,ae) end5、子程序：功能：计算单元的转换矩阵。subroutine cta (co,si,t)dimension t(6,6)call clear(6,6,t)t(1,1)=cot(1,2)=sit(2,1)=-sit(2,2)=cot(3,3)=1.0do 20 i=1,3do 10 j=1,310t(i+3,j+3)=t(i,j)20 continueend7、功能：由整体坐标系下的单刚形成总刚。将局部坐标系单刚所有的的元素置将单刚的上三角中不为0 的元素赋值对单刚的下三角元素赋值将单刚转换成整体坐标系下的单刚Subroutine formk(io,jo,n,ae,a)Dimension a(n,n),ae(6,6)Do 20 I=1,3Ii=io+IDo 10 j=I,3Jj=jo+j10 a(ii,jj)=a(ii,jj)+ae(I,j)20 continuedo 40 I=4,6ii=jo+I-3do 30 j=I,6jj=jo+j-330 a(ii,jj)=a(ii,jj)+ae(I,j)40 continueif (io.lt.jo) thendo 60 I=1,3ii=io+Ido 50 j=4,6jj=jo+j-350 a(ii,jj)=a(ii,jj)+ae(I,j)60 continueelsedo 80 I=4,6ii=jo+I-3do 70 j=1,3jj=io+j70 a(ii,jj)=a(ii,jj)+ae(I,j)80 continueend ifend8、子程序功能：处理结点荷载。SUBROUTINE IOJP(N,P)DIMENSION P(N)!P(N)荷载列阵!NJP 结点荷载数，如果同一个结点在三个方向分别作用有荷载,则应记为三个荷载!JN 荷载作用的结点号;!JD 荷载作用的方向号(1水平力;2垂直力;3弯矩);!PV 荷载值(水平和垂直荷载分别以整体坐标的x和y正向为正;弯矩以逆时针为正); READ(5,*)NJP WRITE(6,(1X,NJP=,I3) NJP IF (NJP.EQ.0) GOTO 200 WRITE(6,(3X,JN,7X,JD,7X,PV) DO 100 I=1,NJPREAD(5,*) JN,JD,PVWRITE(6,(I5,I9,E12.4) JN,JD,PVJ=3*(JN-1)+JDP(J)=P(J)+PV100 CONTINUE200 CONTINUE END9、子程序 功能：处理非结点荷载。 subroutine ftfb(k,n,rl,io,jo,ilt,pv,dx,t,p)common/m2/img(500,2),tf(500,6),ek(500,4)dimension t(6,6),p(n),f(6),f1(6)do 10 I=1,610 f(i)=0.0a=dxb=1.0-a goto (20,40) iltx=dx*rlc第1类荷载类型f(2)=pv*(1.0+2.0*dx)*b*bf(5)=pv-f(2)f(3)=pv*x*b*bf(6)=-pv*a*a*(rl-x)goto 200c 第2类荷载类型f(1)=-pv*b1f(4)=-pv-f(1)goto 200200do 210 j=1,6210 tf(k,j)=tf(k,j)+f(j)do 230 j=1,6f1(j)=0.0do 220 jj=1,6220 f1(j)=f1(j)-t(jj,j)*f(jj)230 continuedo 240 j=1,3ji=io+jjj=jo+jp(ji)=p(ji)+f1(j)240p(jj)=p(jj)+f1(j+3) end根据荷载信息计算固端力将单元的固端力存放在TF（NM，6）数组中按照公式计算整体坐标系下的单元等效荷载列向量。根据单元首尾结点编号IO、JO累加到荷载列阵中的相应位置该子程序的程序框图三总程序PROGRAM ccsIMPLICIT DOUBLE PRECISION (A-H,O-Z ) ,INTEGER (I-N)COMMON /M1/ISE(500,2),XY(500,2)/M2/IMG(500,2),TF(500,6),EK(500,4)DIMENSION IS(100),W(8000),T(6,6),AE(6,6)CHARACTER PN*14,FN*12!NM,NJ,NS,NGE,NME,N 结构参数；!TF(6,6)坐标转换矩阵T;A9(6,6)单元刚度矩阵；DATA W/8000*0/WRITE(*,(A) 输入计算问题名(PN)READ(*,(A) PNCALL FNAME(PN,.DAT,FN)OPEN(5,FILE=FN,STATUS=OLD)CALL FNAME(PN,.OUT,FN)OPEN(6,FILE=FN,STATUS=UNKNOWN)READ(5,*) NM,NJ,NS,NGE,NME!计算各数组对应的一维动态数组的其始位置N=3*NJL1=1 !AI(NGE,2) W(L1)L2=L1+2*NGE !A(N,N) W(L2)L3=L2+N*N !WE(NME,2) W(L3)L4=L3+2*NME !P(N) W(L4)L5=L4+N !SD(NS) W(L5)LL=L5+NSIF(LL.LT.8000) GOTO 10WRITE(*,(A,I5) LL=,LLSTOP COMMON SIZE IS TOO SMALL!10CONTINUE!输入；CALL INPUT(NM,NJ,NS,NGE,NME,IS(1),W(L5),W(L1),W(L3)!处理节点荷载；CALL IOJP(N,W(L4)!处理非结点荷载；READ(5,*) NLM !NLM 非结点荷载个数WRITE(6,(1X,NLM=,I3) NLMIF (NLM.EQ.0) GOTO 200DO I=1,NLMREAD(5,*) K,ILT,PV,DXWRITE(6,205)WRITE(6,210) K,ILT,PV,DXCALL CH1(K,IO,JO,SI,CO,RL)CALL CTA(CO,SI,T)CALL FTFB(k,n,rl,io,jo,ilt,pv,dx,t,W(L4)205FORMAT(2X, K, 7X, ILT, 11X, PV, 15X, DX)210FORMAT(1X, I3, 6X, I3, 6X, E12.4, 8X, F5.3)ENDDO200CONTINUE!形成总刚；DO 100 K=1,NMCALL CH1(K,IO,JO,SI,CO,RL)CALL CTA(CO,SI,T)CALL DKE(K,NME,NGE,RL,W(L1),W(L3)CALL KE(K,T,AE)CALL FORMK(IO,JO,N,AE,W(L2)100CONTINUE!引入支承条件，处理总刚；CALL ASD(NS,N,IS(1),W(L2)!引入支承条件，处理右端项；CALL BSD(NS,N,IS,W(L5),W(L4)!解方程，求结点位移；CALL GAUSS1(N,W(L2),W(L4)!计算杆端力；call cotf(nm,n,w(l4)!输出；CALL OUTPUT(NM,NJ,N,W(L4)END!高斯消去法SUBROUTINE GAUSS1(N,A,B)IMPLICIT DOUBLE PRECISION (A-H,O-Z ) ,INTEGER (I-N)DIMENSION A(N,N),B(N) !系数矩阵消元 DO 40 M=1,N-1 DO 30 I=M+1,N C1=A(M,I)/A(M,M) DO 20 J=I,N A(I,J)=A(I,J)-A(M,J)*C1 20 CONTINUE 30 CONTINUE 40CONTINUE !右端项消元 DO 60 M=1,N-1 C2=B(M)/A(M,M) DO 50 I=M+1,N B(I)=B(I)-A(M,I)*C2 50CONTINUE 60CONTINUE !回代求结点位移 B(N)=B(N)/A(N,N) DO 80 I=N-1,1,-1 DO 70 J=I+1,N B(I)=B(I)-A(I,J)*B(J) 70 CONTINUE B(I)=B(I)/A(I,I) 80CONTINUE END!数据输入SUBROUTINE INPUT(NM,NJ,NS,NGE,NME,IS,SD,AI,WE)IMPLICIT DOUBLE PRECISION (A-H,O-Z ) ,INTEGER (I-N)COMMON /M1/ISE(500,2),XY(500,2)/M2/IMG(500,2),TF(500,6),EK(500,4)DIMENSION IS(NS),AI(NGE,2),WE(NME,2),SXYM(3),SD(NS)DO I=1,NMREAD(5,*) M,(ISE(I,J),J=1,2),(IMG(I,J),J=1,2)ENDDOREAD(5,*) (M,(XY(I,J),J=1,2),I=1,NJ)READ(5,*) (M,(WE(I,J),J=1,2),I=1,NME)READ(5,*) (M,(AI(I,J),J=1,2),I=1,NGE)READ(5,*) NSJJS=0DO 50 JSJ=1,NSJREAD(5,*) ISJ,(SXYM(I),I=1,3)ID=3*(ISJ-1)DO 40 J=1,3IF(SXYM(J).GT.9998.0) GOTO 40JS=JS+1SD(JS)=SXYM(J)IS(JS)=ID+J40CONTINUE50CONTINUEIF(JS.GT.NS) THENWRITE(*,(1X,NS TOO SMALL /1X,NS SHOULD BE,I3) JSSTOPEND IFWRITE(6,100)WRITE(6,110) NM,NJ,NS,NGE