matlab连续梁程序的编制与使用

第三章连续梁程序的编制与使用结构力学中的矩阵位移法是随着电子计算机进入结构力学领域中而产生的一种方法，而Matlab语言正是进行矩阵运算的强大工具，因此，用Matlab语言编写结构力学程序有更大的优越性。本章将详细介绍如何利用Matlab语言编制连续梁结构的计算程序。矩阵位移法的解题思路是将结构离散为单元（杆件），建立单元杆端力与杆端位移之间的关系－单元刚度方程；再将各单元集成为原结构，在满足变形连续条件和平衡条件时，建立整体刚度方程；在边界条件处理完毕后，由整体刚度方程解出节点位移，进而求出结构内力。用矩阵位移法计算连续梁的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，并进行编码：单元编码、结点编码、结点位移编码、选取坐标系。2）建立局部坐标系下的单元刚度矩阵。3）建立整体坐标系下的单元刚度矩阵。4）集成总刚。5）建立整体结构的等效节点荷载和总荷载矩阵6）边界条件处理。7）解方程，求出节点位移。8）求出各单元的杆端内力。实际上，上述步骤也是编制Matlab程序的基本步骤，在求出计算结果后，还可以利用Matlab的绘图功能绘制结构图、内力图、变形图等等。开开始始输输入入初初始始数数据据形形成成原原始始总总刚刚形形成成综综合合结结点点荷荷载载边边界界处处理理解解方方程程求求杆杆端端力力并并输输出出结结束束图图33--11程程序序流流程程图图3.1程序说明%*******************************************************************%矩阵位移法解连续梁主程序%*******************************************************************功能：运用矩阵位移法解连续梁的基本原理编制的计算主程序。基本思想：结点（结点位移）编码默认为从左至右，从1开始顺序进行；杆端弯矩的方向默认为逆时针。荷载类型：可计算结点荷载，每单元作用的跨中集中力和均布荷载。说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。%-----------------------------------------------------------------------------------------------------1程序准备formatshorte%设定输出类型clearall%清除所有已定义变量clc%清屏说明：formatshorte－设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clearall－清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc－清屏，使屏幕在本程序运行开始时%-----------------------------------------------------------------------------------------------------2打开文件FP1=fopen('input.txt','rt');%打开输入数据文件存放初始数据FP2=fopen('output.txt','wt');%打开输出数据文件存放计算结果说明：FP1=fopen('input.txt','rt');－打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来FP2=fopen('output.txt','wt');－打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。%-----------------------------------------------------------------------------------------------------3读入程序控制信息NELEM=fscanf(FP1,'%d',1);%单元个数（单元编码总数）NPOIN=fscanf(FP1,'%d',1);%结点个数（结点编码总数）NVFIX=fscanf(FP1,'%d',1);％约束个数（零位移总数）NFPOIN=fscanf(FP1,'%d',1);％结点荷载个数（作用在结点上集中力偶总数）NFPRES=fscanf(FP1,'%d',1);％非结点荷载数（作用在单元上分布荷载总数）YOUNG=fscanf(FP1,'%f',1);％弹性模量说明：从输入文件FP1中读入单元个数，结点个数，约束个数，结点荷载个数，非结点荷载个数，弹性模量；程序中弹性模量仅输入了一个值，表明本程序仅能求解一种材料构成的结构，如：钢筋混凝土结构、钢结构，不能求解钢筋混凝土－钢组合结构。采用了命令fscanf，其中：’%d’表示读入整数格式；1表示读取1个数。%------------------------------------------------------------------------------------------fprintf(FP2,'\n结构初始数据\n\n');fprintf(FP2,'单元总数=%d结点总数=%d约束总数=%d\n',NELEM,NPOIN,NVFIX);fprintf(FP2,'结点荷载个数=%d非结点荷载个数=%d弹性模量=%12.5g\n',NFPOIN,NFPRES,YOUNG);说明：在输出文件FP2中显示输入的控制信息，便于程序执行完毕后，从输出文件中查找输入错误。采用了命令fprintf，其中：\n为换行标志，表示换一行；%12.5g表示输出12位且有5位小数的实数。括在引号中的信息'单元总数=%d结点总数=%d约束总数=%d\n'为输出到FP2文件中时的格式，其后的变量表NELEM,NPOIN,NVFIX依次将其代表的数值输出到%d所对应的位置。%-----------------------------------------------------------------------------------------------------4调用读取初始数据函数，生成结构信息[LNODS,COORD,FPOIN,FPRES,FIXED]=…ele_INITIALDATA(FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX);%-----------------------------------------------------------------------------------------------------5调用总刚计算函数，生成单刚并集成总刚[HK]=ele_HK(NPOIN,NELEM,LNODS,COORD,YOUNG);%-----------------------------------------------------------------------------------------------------6调用荷载计算函数，由结点荷载与非结点荷载生成总荷载向量列阵[FORCE]=…ele_FORCE(NPOIN,NFPOIN,FPOIN,NFPRES,FPRES,LNODS,COORD);%-----------------------------------------------------------------------------------------------------7调用边界条件处理函数，总刚、总荷载进行边界条件处理[HK,FORCE]=ele_BOUNDARY(NVFIX,FIXED,HK,FORCE);%-----------------------------------------------------------------------------------------------------8方程求解DISP=zeros(NPOIN,1);%建立并初始化结构的结点位移列矩阵DISP=HK\FORCE;%计算出结构所有的结点位移%-----------------------------------------------------------------------------------------------------9调用单元杆端力计算函数，求单元杆端力[FE]=ele_MOMENTS(FP2,NPOIN,DISP,NELEM,LNODS,…COORD,YOUNG,NFPRES,FPRES);%*******************************************************************10关闭输出数据文件fclose(FP2);%*******************************************************************%读取初始数据函数ele_INITIALDATA%*******************************************************************%入口参数：FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX%出口参数：LNODS,COORD,FPOIN,FPRES,FIXED%------------------------------------------------------------------------------------------function[LNODS,COORD,FPOIN,FPRES,FIXED]=…ele_INITIALDATA(FP1,FP2,NELEM,NPOIN,NFPOIN,NFPRES,NVFIX)说明：FP1－文件句柄，指示输入数据文件；FP2－文件句柄，指示输出数据文件；%------------------------------------------------------------------------------------------%读取结构信息LNODS=fscanf(FP1,'%f',[3,NELEM])';说明：建立LNODS矩阵，该矩阵指出了每一单元的连接信息和惯性矩。矩阵的每一行针对每一单元，共计NELEM；每一列相应为单元左结点号（编码）、（编码）、惯性矩。命令中，[3,NELEM]表示读取NELEM行3列数据赋值给LNODS矩阵。显然，LNODS(i,1:3)依次表示i单元的左结点号、右结点号、惯性矩。从这种定义可见，每一单元的惯性矩均为常数，均为等截面直杆。%------------------------------------------------------------------------------------------fprintf(FP2,'单元号左结点号右结点号惯性矩\n');fori=1:NELEMfprintf(FP2,'%3d%3d%3d%10.2e\n',i,LNODS(i,:));end说明：在输出文件FP2中显示输入的单元连接等信息。从1单元到NELEM单元进行循环，并依次输出每一单元的连接等信息。%----------------------------------------------------------------------------------------------------COORD=fscanf(FP1,'%f',[NPOIN])';%坐标：x坐标（共计NPOIN组）说明：建立COORD矩阵，该矩阵用来存储各结点x方向的坐标值。连续梁结构各结点均分布在x轴，以1结点为起始点顺序编码，因此表征各结点位置时仅需存储其x方向的坐标即可。从FP1文件中读取全部结点个数NPOIN的坐标值。COORD(i)表示第i个结