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

- - 82 - 第三章 连续梁程序的编制与使用结构力学中的矩阵位移法是随着电子计算机进入结构力学领域中而产生的一种方法，而 Matlab语言正是进行矩阵运算的强大工具，因此，用Matlab 语言编写结构力学程序有更大的优越性。本章将详细介绍如何利用 Matlab 语言编制连续梁结构的计算程序。矩阵位移法的解题思路是将结构离散为单元（杆件） ，建立单元杆端力与杆端位移之间的关系单元刚度方程；再将各单元集成为原结构，在满足变形连续条件和平衡条件时，建立整体刚度方程；在边界条件处理完毕后，由整体刚度方程解出节点位移，进而求出结构内力。用矩阵位移法计算连续梁的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，并进行编码：单元编码、结点编码、结点位移编码、选取坐标系。2）建立局部坐标系下的单元刚度矩阵。3）建立整体坐标系下的单元刚度矩阵。4）集成总刚。5）建立整体结构的等效节点荷载和总荷载矩阵6）边界条件处理。7）解方程，求出节点位移。8）求出各单元的杆端内力。实际上，上述步骤也是编制 Matlab 程序的基本步骤，在求出计算结果后，还可以利用 Matlab 的绘图功能绘制结构图、内力图、变形图等等。入入 入入入入入入 入入入入入入 入入入入入入入入 入入入入 入入入 入入入入入入入 入入 图图 3-1 程序流程图程序流程图- - 83 - 3.1 程序说明%*% 矩阵位移法解连续梁主程序%* 功能：运用矩阵位移法解连续梁的基本原理编制的计算主程序。 基本思想：结点（结点位移）编码默认为从左至右，从 1 开始顺序进行；杆端弯矩的方向默认为逆时针。 荷载类型：可计算结点荷载，每单元作用的跨中集中力和均布荷载。 说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。%-1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏 说明：format short e 设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all 清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；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); 打开输出数据文件，该文件不存在时，通过- - 84 - 此命令创建新文件，该文件存在时则将原有内容全部删除。该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。%-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 结构初始数据nn);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 调用 读取初始数据 函数，生成结构信息- - 85 - 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=HKFORCE; %计算出结构所有的结点位移%-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) 说明：- - 86 - FP1 文件句柄，指示输入数据文件；FP2 文件句柄，指示输出数据文件；%-% 读取结构信息LNODS=fscanf(FP1,%f,3,NELEM); 说明：建立 LNODS 矩阵，该矩阵指出了每一单元的连接信息和惯性矩。矩阵的每一行针对每一单元，共计 NELEM；每一列相应为单元左结点号（编码） 、 （编码） 、惯性矩。命令中，3,NELEM表示读取 NELEM 行 3 列数据赋值给 LNODS 矩阵。显然，LNODS(i,1:3)依次表示 i 单元的左结点号、右结点号、惯性矩。从这种定义可见，每一单元的惯性矩均为常数，均为等截面直杆。%-fprintf(FP2,单元号 左结点号 右结点号 惯 性 矩 n);for i=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 个结点的 x 坐标。%-fprintf(FP2, 结点定义数据：结点号 x 坐 标 n);for i=1:NPOINfprintf(FP2, %3d %6.2fn,i,COORD(i);- - 87 - end 说明：在输出文件 FP2 中显示输入的结点坐标信息。从 1 结点到 NPOIN 结点进行循环，并依次输出每一结点的坐标信息。%-FPOIN=fscanf(FP1,%f,2,NFPOIN); 说明：建立 FPOIN 矩阵，该矩阵用来存储直接作用在结点上的荷载信息。从 FP1 文件读取 NFPOIN 行 2 列数据，赋值给 FPOIN 矩阵。每一行的第一列表示荷载作用的结点号；第二列表示荷载的数值大小。连续梁结构每一结点仅有 1 个转角位移，相应地直接作用在结点上的荷载为力偶，以顺时针转动为正。FPOIN(i,1)表示第 i 个直接结点荷载作用的结点号，FPOIN(i,2)表示第 i 个直接结点荷载的数值大小。若控制数据 NFPOIN 等于零，则此行命令不执行。%-if NFPOIN0 fprintf(FP2, 结点荷载数据：结点号 M 力偶 n);for i=1:NFPOINfprintf(FP2, %3d %3d %6.2fn,i, FPOIN(i,:);endend 说明：首先判断控制数据 NFPOIN 是否为零，若为零则不需输出数据；若不为零，表明本算例有直接作用在结点上的荷载，需要在输出文件 FP2 中显示输入的结点荷载信息。从第 1 到第 NFPOIN 个直接结点荷载进行循环，并依次输出每个直接结点荷载的顺序号、结点号和荷载大小。%-FPRES=fscanf(FP1,%f,3,NFPRES); 说明：建立 FPRES 矩阵，该矩阵用来存储非结点荷载信息。从 FP1 文件读取 NFPRES 行 3 列数据，赋值给 FPRES 矩阵。每一行的第一列表示非结点荷载作用的单元号；第二列表示荷载类型；第- - 88 - 三列表示荷载的数值大小。连续梁结构的局部坐标系 x 轴以向右为正、y 轴以向上为正、转动方向以顺时针为正。FPRES(i,1)表示第 i 个非结点荷载作用的单元号，FPRES(i,2)表示第 i 个非结点荷载的类型（类型定义在函数 ele_FORCE 中） ，FPRES(i,3)表示第 i 个非结点荷载的数值大小。若控制数据 NFPRES 等于零，则此行命令不执行。%-if NFPRES0 fprintf(FP2, 非结点荷载数据：荷载号 单元号 荷载类型 荷载大小n);for i=1:NFPRESfprintf(FP2, %3d %3d %6.2f %6.2fn,i,FPRES(i,:);endend 说明：首先判断控制数据 NFPRES 是否为零，若为零则不需输出数据；若不为零，表明本算例有非结点荷载，需要在输出文件 FP2 中显示输入的非结点荷载信息。从第 1 到第 NFPRES 个非结点荷载进行循环，并依次输出每个非结点荷载作用的荷载序号、单元号、荷载类型和荷载大小。%-FIXED=fscanf(FP1,%f,NVFIX); 说明：建立 FIXED 矩阵，该矩阵用来存储零位移对应的位移编码信息（即有限制转动的约束使结点的转角位移为零，对应的结点位移编码。对于连续梁结构即固定支座对应的结点编码） 。从 FP1 文件读取 NVFIX 行数据，赋值给 FIXED 列矩阵。每一行表示零位移对应的结点位移编码。FIXED(i)表示第 i 个零位移对应的结点位移编码。若控制数据 NFPRES 等于零，则此行命令不执行。%-if NVFIX0 fprintf(FP2, 零位移约束数据： 位移编号 位移编号n);fprintf(FP2, %3d %3dn,FIXED(:);end- - 89 - 说明：首先判断控制数据 NVFIX 是否为零，若为零则不需输出数据；若不为零，表明本算例有零位移，需要在输出文件 FP2 中显示零位移约束数据。从第 1 到第 NVFIX 个零位移进行循环，依次输出每个零位移对应的结点号。%-return 说明：return 表明主程序终止。%*% 计算总刚矩阵 函数 ele_HK%* % 入口参数：结点数、单元数、单元信息数组、结点坐标、弹性模量% 出口参数：整体刚度矩阵function HK = ele_HK(NPOIN,NELEM,LNODS,COORD,YOUNG)%-HK=zeros(NPOIN,NPOIN); % 生成总刚矩阵并清零% 生成单刚并组成总刚for i=1:NELEM % 对单元个数循环% 调用生成局部单刚(局部坐标) 函数EK=ele_EK(i,LNODS,COORD,YOUNG);HK(i:i+1,i:i+1) = HK(i:i+1,i:i+1)+EK;end 说明：对每一单元顺序循环，依次计算每一单元的局部单元刚度矩阵。EK 为一个临时变量，计算得到 i 单元的单刚 EK 后，通过下一条语句集成入总刚矩阵中。在下一次循环中，则计算 i+1 单元的单刚，并赋值给 EK，再进行集成。显然，随着循环的不断进行，EK 对应着不同单元的单刚，值不断发生改变。单刚到总刚的集成是利用了连续梁结构的特点，以分块的形式放入总刚矩阵中。例如：对于第 i 个单元，其四个单刚元素分别放入总刚矩阵的第（ i,i） 、（i,i+1） 、 （i+1,i） 、 （i+1,i+1）四个位置，因此通过 HK 的下标(i:i+1,i:i+1) 来表示相应的元素位置。由于相邻单元的单刚集成时，有元素叠加的情况，因此相应的总刚矩阵元素的叠加采用了 HK(i:i+1,i:i+1) = HK(i:i+1,i:i+1)+EK 的表达形式。- - 90 - %-return%*% 计算单元刚度矩阵 函数 ele_EK%* % 入口参数：单元号、单元信息、结点坐标、弹性模量% 出口参数：局部单元刚度矩阵 说明：在 ele_HK 函数中调用本函数时，弹性模量的入口参数符号为 YOUNG，而在本函数的输入参数表中，符号为 E。函数调用时可以采用此种用法。%-function EK=ele_EK(i,LNODS,COORD,E)%-NL=LNODS(i,1); %左结点号NR=LNODS(i,2); %右结点号L=COORD(NR)-COORD(NL); %单元长度I=LNODS(i,3); %惯性矩i=E*I/L; 说明：利用已知信息，计算每一单元的线刚度。LNODS(i,1)为 i 单元左结点的结点号赋值给 NL； LNODS(i,2)为 i 单元右结点的结点号赋值给 NR。COORD(NR)为第 NR 个结点的 x 坐标；COORD(NL)为第 NL 个结点的 x坐标；二者的差为单元的长度。LNODS(i,1)为 i 单元的惯性矩。%-% 生成单刚( 局部坐标) 右手坐标系EK = 4*i 2*i; 2*i 4*i;%-return%*- - 91 - % 生成总荷载向量 函数 ele_FORCE%* % 入口参数：结点数，结点荷载个数，结点荷载信息，非结点荷载个数，非结点荷载信息，单元信息，结点坐标% 出口参数：单元固端力左右两端的杆端弯矩function FORCE = ele_FORCE(NPOIN,NFPOIN,FPOIN,NFPRES,FPRES,LNODS,COORD)%-FORCE=zeros(NPOIN,1); % 生成总荷载向量并清零for i=1:NFPOIN % 对结点荷载个数进行循环FORCE(FPOIN(i,1)=FORCE(FPOIN(i,1)+FPOIN(i,2); % 取出结点荷载end 说明：建立并初始化总荷载向量 FORCE，其行数与总刚矩阵、总位移列矩阵相同，即与总结点数 NPOIN 相同（对于连续梁结构） 。对结点荷载个数进行循环，将每个结点荷载的数值放入所在结点号对应的总荷载向量 FORCE 对应的行处FPOIN(i,1)为第 i 个结点荷载所作用的结点号；FPOIN(i,2) 为第 i 个结点荷载的大小；FORCE(FPOIN(i,1)为总荷向量对应结点号处的荷载。%-for i=1:NFPRES % 对非结点荷载个数进行循环F0=ele_FPRES(i,FPRES,LNODS,COORD); % 计算单元固端力% 对单元局部杆端力要进行坐标转换ele=FPRES(i,1); % 取荷载所在的单元号NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号FORCE(NL)=FORCE(NL)-F0(1); %将固端力变成等效结点荷载FORCE(NR)=FORCE(NR)-F0(2); %固端力与等效荷载符号相反end 说明：计算由非结点荷载引起的等效结点荷载，并集成入已有的总荷向量中。对非结点荷载个数进行循环，即有几个非结点荷载循环几次，若没有非结点荷载，则本循环语句不执行。调用计算非结点荷载引起的单元固端力函数，计算出非结点荷载所在单元的两端固端力。- - 92 - FPRES(i,1)为第 i 个非结点荷载所在的单元号赋值给 ele；LNODS(ele,1) 为ele 单元左结点的结点号赋值给 NL；LNODS(ele,2)为 ele 单元右结点的结点号赋值给 NR。FORCE(NL)、FORCE(NR)分别为单元的左、右杆端的杆端力； -F0(1)、-F0(2)为非结点荷载引起的左、右端固端力的负值，即等效结点荷载。FORCE(NL)=FORCE(NL)-F0(1)为直接作用在单元左端结点上的荷载与左端等效结点荷载的叠加；FORCE(NR)=FORCE(NR)-F0(2) 为直接作用在单元右端结点上的荷载与右端等效结点荷载的叠加。%-return%*% 计算非结点荷载引起的单元固端力 函数 ele_FPRES% 符号规定：正方向为：X 向右 Y 向上 M 顺时针%* % 入口参数：荷载序号，非结点荷载信息，单元信息，结点坐标% 出口参数：单元固端力（左右两端固端弯矩）%-function F0=ele_FPRES(i,FPRES,LNODS,COORD) %-ele=FPRES(i,1); % 取荷载所在的单元号G=FPRES(i,3); % 单元荷载大小NL=LNODS(ele,1); NR=LNODS(ele,2); % 单元的左右结点号L=COORD(NR)-COORD(NL); % 单元长度 说明：将已知非结点荷载信息分别赋值给相应的变量。FPRES(i,1)为第 i 个非结点荷载作用的单元号； G=FPRES(i,3)为第 i 个非结点荷载的大小；LNODS(ele,1)为 ele 单元左结点的结点号赋值给 NL；LNODS(ele,2) 为 ele单元右结点的结点号赋值给 NR。COORD(NR)为第 NR 个结点的 x 坐标；COORD(NL)为第 NL 个结点的 x坐标；二者的差为单元的长度。%-F0=0;0; % 单元固端弯矩清零switch FPRES(i,2)- - 93 - case 1 % 横向均布荷载F0(1)=-G*L2/12.0;F0(2)= G*L2/12.0;case 2 % 横向集中力F0(1)=-G*L/8;F0(2)= G*L/8;end 说明：根据荷载类型，求解单元杆端的杆端弯矩。FPRES(i,2)为第 i 个非结点荷载的类型。该类型是在下面的 switch 判断模块中定义的。当 FPRES(i,2)1 时，为横向均布荷载（横向指垂直于杆件轴线） ；当 FPRES(i,2)2 时，为横向集中力。即荷载类型分别为 1、2，这两个数值需要在输入非结点荷载信息时针对荷载的类型进行输入。每种非结点荷载作用下的单元固端力的计算公式是按照结构力学教材得到的。其中，F0(1)表示单元左端的固端弯矩；F0(2)表示单元右端的固端弯矩。%-return%*% 总刚、总荷载进行边界条件处理 函数 ele_BOUNDARY%* % 入口参数：零位移个数，约束信息，总刚矩阵，总荷载矩阵% 出口参数：总刚矩阵，总荷载矩阵%-function HK,FORCE = ele_BOUNDARY(NVFIX,FIXED,HK,FORCE)%-for j=1:NVFIX % 对约束个数进行循环N1=FIXED(j); %找出约束所在的结点位移编码HK(N1,N1)=HK(N1,N1)*1E10; %乘大数法end 说明：采用乘大数法处理边界条件，仅需对总刚元素进行处理即可。FIXED(j)为第 j 个约束所在的结点位移编码（连续梁即结点号）赋值给N1。- - 94 - 对约束个数进行循环，对每 j 个约束处理时，需对总刚矩阵中的第 j 行第 j列的主对角线元素 HK(N1,N1)乘以一个大数，如 。10%-若采用置 0 置 1 法处理边界条件，其相应程序为for j=1:NVFIX N1=FIXED(j); %找出约束所在的结点位移编码HK(:,N1)=0; %将总刚矩阵 HK(:,N1)的第 N1 列所有元素赋值为 0。HK(N1,:)=0; %将总刚矩阵 HK(N1,:)的第 N1 行所有元素赋值为 0。HK(N1,N1)=1; %将总刚矩阵 HK(N1,N1)主对角线元素赋值为 1。FORCE(N1)=0; %将总荷载矩阵 FORCE(N1)的第 N1 行元素赋值为 0end%-return%*% 计算单元杆端力 函数 ele_MOMENTS%* % 入口参数：输出文件句柄，结点数，结构位移列阵，单元数，单元信息，结点坐标，弹性模量，非结点荷载个数，非结点荷载信息% 出口参数：单元杆端力左右两端的杆端弯矩function FE = ele_MOMENTS(FP2,NPOIN,DISP,NELEM, LNODS,COORD, YOUNG, NFPRES, FPRES)%-fprintf(FP2,n 计 算 结 果nn); fprintf(FP2, 输出结构结点位移 结点号 转角n);for i=1:NPOINfprintf(FP2, %3d %11.3e n,i,DISP(i);end 说明：在输出文件 FP2 中显示提示信息：计算结果。通过循环按照结点位移编码（连续梁结构即结点编码）的顺序，在文件FP2 中输出所有结点位移的数值。%- - 95 - EDISP=zeros(2,1); % 单元位移列向量清零fprintf(FP2, 输出单元杆端内力，以顺时针为正n);for i=1:NELEM % 对单元个数进行循环，依次计算每一单元的杆端力for j=1:2N1=LNODS(i,j); % 单元杆端结点号EDISP(j)=DISP(N1); % 取 N1 端的单元位移列向量end 说明：通过循环，依次从结构位移列阵中对号，赋值给第 i 个单元的单元位移向量 EDISP；j=1 时 N1 为左端， j=2 时 N1 为右端。%-EK=ele_EK(i,LNODS,COORD,YOUNG); % 生成局部单刚(局部坐标)FE=EK*EDISP; %计算局部坐标系下的由结点位移产生的杆端力 说明：计算第 i 个单元的单刚 EK；由单刚乘以单元位移得到单元仅由结点位移产生的杆端力。%-for j=1:NFPRESif FPRES(j,1) = i %成立时，当前单元上有非结点荷载F0=ele_FPRES(j,FPRES,LNODS,COORD); %计算单元固端力FE=FE+F0; %考虑由非结点荷载引起的杆端力，直接叠加endend 说明：对非结点个数进行循环。FPRES(j,1)为第 j 个非结点荷载作用的单元号，若等于当前循环的单元号 i，则表明当前单元上有非结点荷载。调用计算单元固端力函数，得到非结点荷载引起的单元固端力。将杆端位移引起的杆端力和非结点荷载引起的单元固端力叠加，得到当前单元的全部杆端力。%-fprintf(FP2, 单元(%d) 杆端力 弯 矩n,i);fprintf(FP2, 左端%3d 结点 %7.2fn,LNODS(i,1),FE(1);fprintf(FP2, 右端%3d 结点 %7.2fn,LNODS(i,2),FE(2); 说明：输出单元的结点号，对应的杆端弯矩。- - 96 - end %对单元循环结束。%-return3.2 程序应用举例【例 1】利用连续梁程序计算图 2-9 所示的连续梁结构。 1EI图 2-9%-输入数据文件 input.txt 为：3 4 1 1 2 1.01 2 1.02 3 1.03 4 1.00.06.012.017.02 6.01 2 -35.03 1 -6.04EImkN6II353 5m63kN6- - 97 - %-说明：第一行：读入程序控制信息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); 弹性模量第二、三四行：读入单元连接信息：LNODS=fscanf(FP1,%f,3,NELEM); %单元号，左结点号，右结点号，惯性矩第五、六、七、八行：读入结点坐标COORD=fscanf(FP1,%f,NPOIN); %结点坐标值第九行：读入结点荷载信息FPOIN=fscanf(FP1,%f,2,NFPOIN); %结点号，结点荷载数值（以顺时针为正）第十、十一行：读入非结点荷载信息FPRES=fscanf(FP1,%f,3,NFPRES); %荷载号，单元号，荷载类型，荷载大小%荷载类型：1 为均布荷载，2 为集中力第十二行：读入零位移信息FIXED=fscanf(FP1,%f,NVFIX); %零位移约束所在的结点位移编码%-输出数据文件 output.txt结构初始数据- - 98 - 单元总数=3 结点总数=4 约束总数=1 结点荷载个数=1 非结点荷载个数=2 弹性模量= 1 单元定义数据：单元号 左结点号 右结点号 惯 性 矩 1 1 2 1.00e+000 2 2 3 1.00e+000 3 3 4 1.00e+000 结点定义数据：结点号 x 坐 标 1 0.002 6.003 12.004 17.00结点荷载数据：结点号 M 力偶 1 2 6.00非结点荷载数据：荷载号 单元号 荷载类型 荷载大小1 1 2.00 -35.002 3 1.00 -6.00零位移约束数据： 位移编号 位移编号4 计 算 结 果输出结构结点位移 结点号 转角1 5.597e+001 2 -3.320e+001 3 1.607e+001 4 -2.366e-009 输出单元杆端弯矩 （顺时针为正）单元(1)杆端力 弯 矩左端 1 结点 0.00右端 2 结点 22.78单元(2)杆端力 弯 矩左端 2 结点 -16.78右端 3 结点 -0.35单元(3)杆端力 弯 矩左端 3 结点 0.35- - 99 - 右端 4 结点 18.93%-根据单元杆端力计算结果，可得结构的弯矩图为：3.4 程序的改进要点上述连续梁程序反映了矩阵位移法的基本思路，可以计算简单的连续梁问题。在熟练掌握了程序的编制与使用后，可在以下几方面对程序进行改进，以加深对矩阵位移法及 MATLAB 语言编程的理解：1、本程序的单元信息输入、结点坐标输入均分别隐含要求从 1 单元、1 结点按编码的序号顺序输入，考虑如何改进使相关信息的输入可以更加灵活。2、本程序的弹性模量仅能输入一个数值，意味着程序仅能计算由同种材料构成的结构。考虑如何改进使程序可计算由不同材料构成的组合结构。3、本程序仅能计算两种非结点荷载类型，考虑如何编制跨中作用集中力偶、单元任意位置作用横向集中力、单元部分范围内作用均布荷载等类型的程序。4、本程序的单元杆端力仅能计算杆端弯矩，考虑如何编制计算杆端剪力的程序5、连续梁结构由于局部坐标系与整体坐标系相同，因此未涉及到坐标转换的问题，考虑如何编制坐标转换的程序。6、连续梁结构由于每个结点仅有一个转角位移，所以程序中单刚、总刚、荷载矩阵等多处均采用结点号来替代结点位移号，进行下表索引和循环控制。但对于其他结构如平面刚架，每个结点一般有三个结点位移，则程序的编制过程中需要严格分清结点号和结点位移号。考虑如何编制平面刚架的计算程序。综上所述，本章的连续梁程序体现了如何将矩阵位移法的计算方法和过程用 MATLAB 程序语言表达出来，重点放在