安徽建筑工程学院计算结构力学8课件.ppt

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