1、计算结构力学,第八章 结构分析程序设计,8-1 概述:程序设计的基本概念与要点,至此,我们已完成了结构矩阵分析杆系有限元法基本原理的教学内容,本章主要介绍计算机实现过程。程序设计:当人们需要计算机完成科学计算,数据处理等计算工作时,必须事先恰当地安排好计算步骤,每一步的计算内容以及在什么条件下进行这一部分计算,这就是所谓编制计算机程序。,程序设计应注意以下几点： 保证程序的正确性，可通过考题校验 使程序具有高效率,并尽可能降低费用,求解方程组最费机时(80%左右),应设计再解功能 使程序便于调试、修改、扩充和完善,既要有通用性,又要留下可扩充修改的窗口,从结构矩阵分析原理到最终计算机实现解决具
2、体问题,主要有以下过程模块,用框图可表示为：,8-2 结构分析程序设计的框图设计,根据结构方程组的建立与求解来进行程序设计 结构方程组的建立与求解是结构分析的核心问题,如建立结构方程组的有限元方法,解线性代数方程组的消元分解法,这两个方法一经确定,程序设计的大致思路与过程也就基本确定了。,1. 编程要点,由单元定位向量组织整体流程图的运行实施从形成结构方程组K=P过程中K与P的形成,到计算结构内力和反力的过程,都离不开MW的组织。 应具有良好的通用性本程序的编制特点是利用特殊结点的约束信息，可模拟六种类型的杆系结构，故具有广泛的适用性。,2. 设计流程图,3. 框 图 设 计,4. 本程序设计
3、的模块功能介绍,模块：具有单一的独立的功能块,由子程序或自定义函数所组成。,模块可根据其功能进一步划分,依次分解成较低级的模块,模块之间通过调用而组成一个协同的程序；这种可通过自上而下进行分解,并可通过自上而下的调用,一级一级地组成程序是程序设计的重要方法。,各个模块的功能在很大程度上是独立的,因而不同的模块可以由不同的人来完成。例如,不太了解矩阵力学的人也可以设计消元分解及正消回代的子程序SUB. NXFJ。 模块的相互独立性不仅方便程序设计,也方便程序调试。调试时自下而上一块一块地进行。这时主要着重调试模块接口和上一级模块,而不必调试已通过的模块。 亦可以根据模块功能组成其它程序。,本程序
4、的模块设有三个级别:0、1、2 0级表示各个程序均可调用 1级仅供2级调用。,5. 本程序的静力计算功能 (1)结构形式:可对六种类型的杆系进行计算分析 连续梁桁架排架刚架框排架梁、桁组合结构。可由JTX(4,NJT)进行分类。 (2)材料:各向同性,按EI(EA)进行分组,分组数为NAI,(3)荷载类型 结点荷载信息 NPJ:受载结点数 需输入结点荷载信息数组PJZ(NPJ,2),NPJ行,2列: 第一列输JD.x,第二列输大小(与坐标一致为正，反之为负)。 JD是结点号，x是方向(1为X向, 2为Y向, 3为z向)。,单元荷载信息 单元荷载信息数组PMZ(NPM,3): NPM:受载单元数
5、 NPM行,3列,第一列输M.L:第M单元，第L类荷载； 第二列输荷载位置(距始端距离x) 第三列输荷载大小(与坐标一致为正，反之为负),L类：共六类,见讲义附表(P63)或参考教材表 5-1(P131),程序PSTDY的子程序SUB.DXL中留有 用户入口,可接入这六种以外的荷载。,8-3 单元定位向量的主线作用,从程序设计框图可以看出:程序设计的每一个环节都离不开单元定位向量,故它在程序设计中起到组织者的作用。我们称之为主线作用。,回顾:单元定位向量是按单元结点编号顺序由结点的未知量编号所组成的一个列向量。可由JW(3,NJ)直接生成,其作用主要有：,1.装备结构刚度矩阵 （1）按方阵存贮
6、SUB.KJX1DO 10 I=1,6L=MW(I)IF (L.LE.0) GOTO 10DO 20 J=1,6K=MW(J)IF (K.LE.0) GOTO 20ZK(L,K)=ZK(L,K)+DK(I,J) 20 CONTINUE 10 CONTINUE,思考题：若形成上三角阵如何改动 （2）一维变带宽上三角按行存贮SUB.KJXDO 10 I=1,6 L=MW(I) IF (L.LE.0) GOTO 10 II=KD(L) DO 20 J=1,6 K=MW(J) IF (K.LE.L) GOTO 20 IJ=II-L+K ZK(IJ)=ZK(IJ)+DK(I,J) 20 CONTINUE
7、 10 CONTINUE,2.形成P 见PSTDY中SUB.YDX 3.单元位移的形成 见PSTDY中SUB.YWY 由(存在P中，存在D(NE,6)中 4.FjRj亦由MW的第j个分量是否为零来判断,见PSTDY中SUB.QDL,5.在KD数组中求带宽亦用到MW 见PSTDY中SUB.QKD 带宽公式：NDK=单元两端未知量编号最大差值+1 6.更为重要的是单元定位向量还体现了单元间的相互联接，以及对结构边界条件的处理，如主从关系、无效未知量的处理等，可通过对特殊结点的约束信息数组JTX(4,NJT)来模拟实际结构中的复杂关系。换句话说，我们在程序的功能中所提到的六种类型的杆系结构，程序设计
8、中最后的区别形式就是MW，这点尤为重要。,对各种类型的结构，我们设计了统一的单元刚度矩阵形式，由MW来直接装配总刚！ 六种类型杆系结构的JTX数组 1.连续梁,2.桁架,3.排架,4.框架,5.框排架,6.框排架,根据各程序模块的功能，本程序设计了下列控制变量与循环变量：NE单元总数NJ结点总数NJT特殊结点数NJZ支座结点数NAIEA或EI的分组数NPJ结点荷载数NPM单元荷载数,8-4 变量与数组设计,N未知量总数M单元序号LOAD荷载分组序号I、J单元刚度矩阵的行列号L，K结构刚度矩阵的行列号II，LL对角线元素地址N1-KD数组的个数，NZY结构刚度矩阵的元素数,整型数组有JH(2，N
9、E)单元两端的结点号MW(6)单元定位向量JW(3，NJ)结构结点未知量编号JTX(4，NJT)特殊结点信息，JZH(NJZ)支座维点号NLX(2，NL)每组荷载信息，即NPJ数和NPM数NAI(1，NE)截面特性分组号KD(NI)结构刚阵一维存贮时主元素地址；,双精度实型数组X(NJ)结点坐标值Y(NJ)结点坐标值SL(NE)单元长度CX(NE)单元的cosSY(NE)单元的sinEA(NAI)单元的EAEI(NAI)单元的EIXS(7)整体坐标系下各单元的7个常数，XSA(NE，7)整体坐标系下各单元的7个常数,DK(6，6)单元刚度矩阵，T(6,6)坐标变换矩阵，ZK(NZY)结构刚度矩
10、阵，一维存贮P(N)荷载列阵D(NE，6)单元两端的结点位移，后存结点力F(6)单元结点力FE(6)等效结点力FG(6)整体坐标系下的单元结点力DG(6)整体坐标系下单元结点位移FLZ(NJZ，3)支座反力PMZ(NPM，3)单元荷载信息PJZ(NPJ，2)结点荷载信息,在采用FORTRAN语言编制程序时，一定要摘清楚数据如何传递。鉴于FORTRAN语言的模块化性质，各程序的数据一般可通过下列三种方式进行传递： 哑实结合 COMMON块 数据文件 方法还可运用在机器设备与外部设备(终端)的数据传递，如本程序设计就采用2个OPEN语句，建立了输入、输出数据文件。 我们在这里介绍的程序设计各子程序
11、间的数据传递均采用方法 。,8-5 数据传递与动态数组设计,子程序的一般形式为 SUBROUTINE QJW(NJ，NJT，JTX，JW，N) RETURN END 其中括号内的NJ，N即为形式参数或称为虚拟变量，或称哑元。它可以是变量名字、数组名字或数组元素。,1、哑实结合的数据传递方式,子程序中的形式参数没有确切的数值，这就是虚拟变量的由来。只有在调用该子程序时，才对形式参数赋值，或赋予其实在的存贮空间。 如：CALL QJW(NJ，NJT，JTX，JW，N) 这时NJ，N称为实在参数(实元)。注意：形式参数和实在参数的类型应一致，个数应相等，但参数名可以不相同。,形式参数的作用可分为两种
12、： 一种是从主程序或其它子程序来赋值的，也就是通过哑实结合接收从外面输入的数据，作为本子程序计算的依据，好象是加工厂的原料一样 另一种是本程序模块计算的结果，通过哑实结合传递到调用处，这是向外传递的数据，好象是加工厂的产品一样。如： 子程序：SUB. QJW(NJ，NJT，JTX，JW，N) 主程序：CALL QJW(NJ，NJT，JTX，JW，N),动态数组又称为可调数组，如：JH(2,NE)、JW(3,NJ)JTX(4,NJT)、ZK(NZY)等，由于这里NE、NJ、NJT，NZY等均为数组变量，没有确切的数字，即数组的大小未能得到确切的定义，因而可调数组在主程序中是不允许出现的。 在主程
13、序中只能出现确切定义的数组，如我们在前二个大作业里所介绍的数组如JH(2,20)、ZK(50,50)等，这样才能在在DIMENSION语句中予以确切定义。,2、主程序中动态数组的设计,根据子程序的形式参数(哑元)的定义，在子程序中可以出现动态数组，但必须在调用时进行哑实结合，才能进行运算。 由于在实际计算中大多数数组的大小是随具体问题的不同而变化的，但在主程序中又不允许出现动态数组，这就给我们在DIMENSION语句中如何定义数组带来不少困难：既受到计算机内存的限制又应使数组有足够的存贮空间，这对于在微机上解决大型工程问题显得更重要。,我们知道，二维或高维数组在计算机内部都是按列存放的，即在计
14、算机内部都是按一维数组 的方式来存贮的，这表示数组变量之间存在一定的关系，如对于：JH220数组： JH(2，3)JH(6)，JH(1，4)JH(7)， 这说明数组一经定义，这种关系便确定，计算机立即“了解”。此外，一个数组变量(下标变量)的下标值加1，就是紧跟在它后面的下标值，这叫做数组变量(下标变量)的后继函数，计算机的处理功能使程序会自动按后继函数找到下一个元素，这叫做下标的自动后继性质。,如在主程序DIMENSION语句中定义说明了JH(2，20)，则在调用时实元用JH或JH(1，1)或JH(1)调用的效果是相同的，如在主程序中用一维定义说明了JH(40)，则用JH或JH(1)调用效果
15、亦是相同。至于调用的JH具体形式，则由子程序对可调数组JH(2，NE)进行说明确定，换句话说，利用下标的自动后继性质，采用首元素调用，即可在主程序中出现动态数组。 根据算例的需要，我们在主程序中仅开设了整、实二个大数组：整：IA(1000)，暂定其大小NIA=1000实：A(10000)，暂定其大小NA=10000,然后按照程序中所要出现数组名字顺序，定出每个数组的第个元素地址，由所输入的变量或巳确定的变量确定这一点并不难，参考下图： 整：,实：,JH(2，NE)的第一个元素地址为IA(1) NLX(2，NL)的第一个元素地址为IA(K1)，这里 (K1)=1+2*NE JW(3，NJ)的第一
16、个元素地址为IA(K2)，这里 (K2)=K1+2*NL 对于各双精度实型数组在A中第一个元素地址亦可由图方便推出。如将某数组的第一个元素作为实元调用，通过哑实结合过程中的下标自动后继，当子程序哑元表中相应的哑元为动态数组时，便可完全按照子程序中该数组的DIMENSION语句的定义在主程序相应数组中得到反映，即在主程序中实现了动态数组。,由此可知，只要找出控制各数组大小的一些变量，即可确定各数组的第一个元素。仔细研究这些变量，发现其中某些变量(如N，N1，NZY)等可通过另一些变量由程序计算确定，我们将后者称为主控变量，须按其出现的次序在程序中首先输入。本程序的主控变量为：NE、NJ、NJT、
17、NJZ、NL、NAI。为促进同学们今后的工程应用能力，我们这里所提供的程序PSTDY并有动力分析内容，故主控变是还增加一个MJ。,MJ：拟求振型数。在PSTDY中，我们将MJ采用屏幕输入，并兼作静、动力分析的程序运行控制，在屏幕的提示下，若键人MJ=0仅作静力计算MJ=3做完静力计算后，再求前三阶频率和振型，具体请参阅算例1，不另赘述。另外，在主程序中还备有这两大类型数组实际使用的元素个数显示，如发生IA或A的溢出现象，只需在内存允许的范围内修改这四条语句，见主程序的001、002、004、005句。,1、根据设计框图，各子程序及其功能汇总如下： 22 KJX：形成结构刚阵ZK(NZY) 23
18、、25 NXFJ：消元分解法解线代方程组23：刚阵ZK的消元25：对P的正消回代 24 YDX：形成P 26 YWY：打印杆端位移 27 QDL：计算F及R,8-6 源程序设计,28(补充) MJX：形成结构质量的矩阵 29(补充) DYNA：逆迭代法计算结构前MJ阶振型和频率 11 QJW：结点未知量编号数组 12 QJZH：形成支座结点号 13 DCH：形成单元常数 14 QKD：形成主元地址数组 15 QXS：形成单刚系数数组 16 DKX：形成单刚矩阵 17 PGP：叠加荷载，形成P的具体计算，,18 DXL：形成FE 01 ZER01：对向量充零 02 ZER02：对二维数组充零 0
19、3 JZZ：矩阵转置 04 QMW：形成单元定位向量。 这里的顺序由二位数表示，十位数上有0，1，2级，其中高级别的可调用低级别的模块，个位数则表示被调用的顺序。,由 OPEN(1，FILE=QAZ.TXT)建立数据文件 READ(1，*) NE，NJ，NJT,NJZ，NL，NAI(见主程序语句标号003) READ(1，*) JH，NLX，JMH，JTX(见主程序语句标号006) READ(1，*) EA，EI，X，Y(见程序语句标号007) READ(1，*) PJZ(NPJ，2)(见SUBYDX语句标号008) READ(1，*) PMZ(NPM，3)(见SUBYDX浯句标号009),2
20、、数据文件的形成,在PSTDY程序中，我们用OPEN语句建立了数据输入文件QAETXT及结果输出文件 FCADTXT。全部的原始数据，除拟求振型数MJ兼作静、动力计算的运行控制从屏幕输入外，均可在QAETXT文件上从终端读人。 全部的数据均采用自由格式，由READ(1，*)语句读人，这里“1”表示通道号，“*”表示自由格式。自由格式要求两个数之间用“，分隔，且每个数所占的位数不限，但要与对应变量的隐含数型一致。,3. 算例,1数据文件可在每行开始直接输入数据，不用说明符号“C”，也没有语句标号区及续行区，但每行不应超过72列 2每个输入语句的结尾不要加任何标点符号，但实型数据的小数点除外 3下
21、一个语句不要与上一个输入语句接着输入，而应别起一行,建立数据文件时还要注意以下几点,算例1作图示刚架的静力计算，并求出其基本频率与振型。,任取,l=4m,m=3.0102kg/m P=2KN,E=2.1106KN/m2 I=4105cm4=410-3m4EI=8.4103KNm2,则,各结点标号与单元划分如图示。,1第一个输入语句在主程序003句，要求输入控制变量6个：NE， NJ， NJT， NJZ， NL， NAI3， 4， 4， 2， 1， 2 注意，这六个整形数成为一行，句尾不加任何符号。,2第二个输入语句在主程序006句，要求输入整型数组IA，实际上输入JH，NLX，JMH，JTX，
22、其中JH有2*NE=6个数，NLX有2*NL=2个数，JMHE有NE=3个数，JTX有4*NJT=16个数，总共29个整型数，可分为两行，但这两行之间要用逗号“，”隔开。 3,1,4,2,1,2,1,0,1,1,2,1,0,10001,0,2,1001,10001,0,3,1,1,1,4,1,1,1 注意：这里特殊节点约束信息JTX(4，NJT)的填写是一项非常细致的工作，需根据具体力学模型，参照17中的有关说明，认真填写。,3第3个输入语句在主程序007句，要求输入双精度型数组A，实际上输入EA，EI，X，y。亦采用自由格式，EA有2个数，EI也有2个数，X有4个数，Y也有4个数，总共有12
23、个数，可一行输入： 0,0,8400,33600,0,8,0,8,4,4,0,0 4在主程序对荷载的循环中，要输入每一组荷载的数据，PSTDY程序仍采用通常的荷载信息集约方法，在集中荷载作用时，要给出节点荷载的个数，荷载作用在哪号节点的哪个位移方向以及荷载的大小，这些信息可以定义一个数组PJZ(NPJ，2)来存放。,PJZ(I，1)可填实型数JD.x 其中：JD一荷载作用的节点号；x一荷载的作用方向，可在总体坐标系下考虑。x=1，荷载作用沿X方向；x=2，荷载作用沿Y方向,x=3，绕Z轴正向作用的力矩。PJZ(1，2)填荷载的大小，与坐标方向一致者为正，I表示节点荷载序号。在单元荷载作用下，P
24、STDY的静力计算可解决六种类型的单元荷载问题，我们只要给出如下信息：,(1)在哪个单元M上有什么类型L的荷载作用； (2)在右手系下荷载离单元左节点的距离； (3)荷载大小Q，其符号规定与坐标方向一致为正。 于是也可定义一个数组PMZ(NPM，3)来存放上述信息。PMZ(I，1)填M.L；PMZ(I，2)填x；PMZ(I，3)填Q；,其中， L=1，表示左端有局部均布荷载Q作用； L=2，表示离左端x处有集中荷载Q作用； L=3，表示离左端x处有集中力矩Q作用； L=4，表示左端有局部三角形荷载作用，其x处为 Q L=5，表示左端有均布轴力荷载Q作用； L=6，表示离左端x处有集中轴力Q作用
25、；,在子程序YDX中，有两个输入语句，应分别输入PJZ与PMZ。本例只有一组结点荷载NL=1，没有单元荷载NPM=0，结点荷载数NPJ=1，所以第四个输入语句在子程序SUBYDX中标号为008句，输入结点荷载信息PJZ，共2*NPJ=2个数，也占行：1.1，2000 注意，由于NPM=0，单元荷载输入信息PME(NPM，3)的输入语句009，程序自动跳转，该语句可以不输入。,5至此，静力计算的所有信息输入完毕，如不做动力计算，可直接从屏幕上输“0”，程序运行结束，并输出结果。如需要进行动力计算时，应输入所求振型数MJ，从屏幕上输入。本例MJ=1，可直接输入“1”。当MJ0时，还应在QAZTXT
26、文件上继续输入动力计算的有关数据，这样，第五个输入语句在子程序SUBMJX的010句，输入ERM(NAI)杆件质量线密度数据，本题有NAI=2个数，也占一行：300，450,解 根据图55，按程序中输入语句的顺序依次填写数据。由算例，可知： 1第一个输入语句标号003，要求输人控制变量6个：NE，NJ，NJT，NJZ，NL，NAI9， 8， 2， 2， 1， 2 这6个数恰好一行，在文件中叫做一个记录。,算例2 做54节例3三层刚架的计算。,2第二个输入语句标号006，要求输入整型数组IA，实际输入JH，NLX，JMH，JTX。这时也按自由格式输入，JH有2*NE=18个数，NLX有2*NL=
27、2个数，JMH有NE=9个数，JTX有4*NJT=8个数，总共37个数，分为3行 3，1，4，2，5，3，6，4，7，5，8，6，1，2，3，4，5，6，0，3，1，1，1，1，1，1，2，2，2，7，1，1，1，8，1，1，1 注意上一个不要与下一个输入语句接着输入，而应另起一行。,3第三个输入语句标号为007，要求输入双精度型数组A，实际上输入EA，EI，X，Y。这时亦采取自由格式，EA有2个数，EI也有2个数，X有8个数，Y也有8个数。总共20个数，占2行： 468000，390000，14040，14040，8128，0，6，0，6，0，6，0，6，10，10，7，7，4，4，0，0,
28、4在对荷载循环中，要输入每一组荷载的数据。在子程序YDX中有两个输入浯句。本例NL=1，没有结点荷载NPJ=0，只有单元荷载NPM=3，所以第四个输入语句标号为009，输入单元荷载信息PMZ，共3*NPM=9个数，这个输入语句只占一行： 7.1,6,-3,8.1,6,-3,9.1,6,-3 这样建立起来的数据文件为“QAZ.TXT“。 部分参考答案：支座反力,算例3作图示连续梁的静力计算,解：各结点标号及单元划分如图示。依前所述，本题数据文件QAE.TXT为： 5,6,6,5,1,1 1,2,2,3,3,4,4,5,5,6,1,3, 1,1,1,1,1,1,1,1,1,2,10001,1,0,3,10001,1,0,4,10001,1,0,5,10001,1,0,6,10001,0,0 0,24,0,4,10,16,20,22,0,0,0,0,0,0 6.2,-20 2.2,2,-90,3.1,6,-20,4.2,2,-60 给出部分结果，并与文5解作比较(设使梁上侧受拉的弯矩为正)。 连续梁部分截面弯矩,
