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

第三章连续梁程序的编制与使用

入结构力学领域中而产生的一种方法，而Matlab

语言正是进行矩阵运算的强大工具，因此，用

Matlab语言编写结构力学程序有更大的优越性。本

章将详细介绍如何利用Matlab语言编制连续梁结

构的计算程序。

矩阵位移法的解题思路是将结构离散为单元

（杆件），建立单元杆端力与杆端位移之间的关系

－单元刚度方程；再将各单元集成为原结构，在满

足变形连续条件和平衡条件时，建立整体刚度方

程；在边界条件处理完毕后，由整体刚度方程解出

节点位移，进而求出结构内力。

用矩阵位移法计算连续梁的步骤如下：

1）整理原始数据，如材料性质、荷载条件、约

束条件等，并进行编码：单元编码、结点编

码、结点位移编码、选取坐标系。

2）建立局部坐标系下的单元刚度矩阵。

3）建立整体坐标系下的单元刚度矩阵。

4）集成总刚。

5）建立整体结构的等效节点荷载和总荷载矩

阵

6）边界条件处理。

7）解方程，求出节点位移。

8）求出各单元的杆端内力。

实际上，上述步骤也是编制Matlab程序的基本

步骤，在求出计算结果后，还可以利用Matlab的绘图功能绘制结构图、内力图、

变形图等等。

图3-1程序流程图

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'); －打开输出数据文件，该文件不存在时，通过此

命令创建新文件，该文件存在时则将原有内容全部删除。该文件设置为可写格式，可在程序执行过程中向输出文件写入数据。

%----------------------------------------------------------------------------------------------------- 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');

for i=1:NELEM

fprintf(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:NPOIN

fprintf(FP2,' %3d %6.2f\n',i,COORD(i));

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 NFPOIN>0

fprintf(FP2,' 结点荷载数据：结点号M力偶\n');

for i=1:NFPOIN

fprintf(FP2,' %3d %3d %6.2f\n',i, FPOIN(i,:));

end

end

●说明：

首先判断控制数据NFPOIN是否为零，若为零则不需输出数据；若不为零，表明本算例有直接作用在结点上的荷载，需要在输出文件FP2中显示输入的结点荷载信息。

从第1到第NFPOIN个直接结点荷载进行循环，并依次输出每个直接结点荷载的顺序号、结点号和荷载大小。

%----------------------------------------------------------------------------------------------------- FPRES=fscanf(FP1,'%f',[3,NFPRES])';

●说明：

建立FPRES矩阵，该矩阵用来存储非结点荷载信息。

从FP1文件读取NFPRES行3列数据，赋值给FPRES矩阵。

每一行的第一列表示非结点荷载作用的单元号；第二列表示荷载类型；第三列表示荷载的数值大小。连续梁结构的局部坐标系x轴以向右为正、y轴以向上为正、转动方向以顺时针为正。

FPRES(i,1)表示第i个非结点荷载作用的单元号，FPRES(i,2)表示第i个非结点荷载的类型（类型定义在函数ele_FORCE中），FPRES(i,3)表示第i个非结点

荷载的数值大小。

若控制数据NFPRES等于零，则此行命令不执行。

%------------------------------------------------------------------------------------------

if NFPRES>0

fprintf(FP2,' 非结点荷载数据：荷载号单元号荷载类型荷载大小\n');

for i=1:NFPRES

fprintf(FP2,' %3d %3d %6.2f %6.2f\n',i,FPRES(i,:));

end

end

●说明：

首先判断控制数据NFPRES是否为零，若为零则不需输出数据；若不为零，表明本算例有非结点荷载，需要在输出文件FP2中显示输入的非结点荷载信息。

从第1到第NFPRES个非结点荷载进行循环，并依次输出每个非结点荷载作用的荷载序号、单元号、荷载类型和荷载大小。

%----------------------------------------------------------------------------------------------------- FIXED=fscanf(FP1,'%f',NVFIX)';

●说明：

建立FIXED矩阵，该矩阵用来存储零位移对应的位移编码信息（即有限制转动的约束使结点的转角位移为零，对应的结点位移编码。对于连续梁结构即固定支座对应的结点编码）。

从FP1文件读取NVFIX行数据，赋值给FIXED列矩阵。

每一行表示零位移对应的结点位移编码。

FIXED(i)表示第i个零位移对应的结点位移编码。

若控制数据NFPRES等于零，则此行命令不执行。

%----------------------------------------------------------------------------------------------- if NVFIX>0

fprintf(FP2,' 零位移约束数据：位移编号位移编号\n');

fprintf(FP2,' %3d %3d\n',FIXED(:));

end

●说明：

首先判断控制数据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的表达形式。

%----------------------------------------------------------------------------------------------------- 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

%******************************************************************* % 生成总荷载向量函数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

●说明：

计算由非结点荷载引起的等效结点荷载，并集成入已有的总荷向量中。

对非结点荷载个数进行循环，即有几个非结点荷载循环几次，若没有非结点荷载，则本循环语句不执行。

调用计算非结点荷载引起的单元固端力函数，计算出非结点荷载所在单元的两端固端力。

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)

case 1 % 横向均布荷载

F0(1)=-G*L^2/12.0;

F0(2)= G*L^2/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。

对约束个数进行循环，对每j个约束处理时，需对总刚矩阵中的第j行第j

10。

列的主对角线元素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行元素赋值为0 end

%---------------------------------------------------------------------------------------------------- return

%******************************************************************* % 计算单元杆端力函数ele_MOMENTS

%******************************************************************* % 入口参数：输出文件句柄，结点数，结构位移列阵，单元数，单元信息，结点坐标，弹性模量，非结点荷载个数，非结点荷载信息

% 出口参数：单元杆端力—左右两端的杆端弯矩

function [FE] = ele_MOMENTS(FP2,NPOIN,DISP,NELEM, LNODS,…

COORD, YOUNG, NFPRES, FPRES)

%---------------------------------------------------------------------------------------------------- fprintf(FP2,'\n 计算结果\n\n');

fprintf(FP2,' 输出结构结点位移结点号转角\n');

for i=1:NPOIN

fprintf(FP2,' %3d %11.3e \n',i,DISP(i));

end

●说明：

在输出文件FP2中显示提示信息：计算结果。

通过循环按照结点位移编码（连续梁结构即结点编码）的顺序，在文件FP2中输出所有结点位移的数值。

%---------------------------------------------------------------------------------------------------- EDISP=zeros(2,1); % 单元位移列向量清零

fprintf(FP2,' 输出单元杆端内力，以顺时针为正\n');

for i=1:NELEM % 对单元个数进行循环，依次计算每一单元的杆端力for j=1:2

N1=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:NFPRES

if FPRES(j,1) == i %成立时，当前单元上有非结点荷载

F0=ele_FPRES(j,FPRES,LNODS,COORD); %计算单元固端力

FE=FE+F0; %考虑由非结点荷载引起的杆端力，直接叠加end

end

●说明：对非结点个数进行循环。

FPRES(j,1)为第j个非结点荷载作用的单元号，若等于当前循环的单元号i，则表明当前单元上有非结点荷载。

调用计算单元固端力函数，得到非结点荷载引起的单元固端力。

将杆端位移引起的杆端力和非结点荷载引起的单元固端力叠加，得到当前单元的全部杆端力。

%------------------------------------------------------------------------------------------

fprintf(FP2,' 单元(%d)杆端力弯矩\n',i);

fprintf(FP2,' 左端%3d结点%7.2f\n',LNODS(i,1),FE(1));

fprintf(FP2,' 右端%3d结点%7.2f\n',LNODS(i,2),FE(2));

●说明：输出单元的结点号，对应的杆端弯矩。

end %对单元循环结束。

%---------------------------------------------------------------------------------------------------- return

3.2 程序应用举例

【例1】利用连续梁程序计算图2-9所示的连续梁结构。1

EI

= Array

图2-9

%----------------------------------------------------------------------------------------------------- 输入数据文件input.txt为：

3 4 1 1 2 1.0

1 2 1.0

2 3 1.0

3 4 1.0

0.0

6.0

12.0

17.0

2 6.0

1 2 -35.0

3 1 -6.0

4

%----------------------------------------------------------------------------------------------------- 说明：

第一行：读入程序控制信息

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

结构初始数据

单元总数=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.00

2 6.00

3 12.00

4 17.00

结点荷载数据：结点号M力偶

1 2 6.00

非结点荷载数据：荷载号单元号荷载类型荷载大小

1 1 2.00 -35.00

2 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

右端4结点18.93

%----------------------------------------------------------------------------------------------------- 根据单元杆端力计算结果，可得结构的弯矩图为：

22.78

3.4 程序的改进要点

上述连续梁程序反映了矩阵位移法的基本思路，可以计算简单的连续梁问题。在熟练掌握了程序的编制与使用后，可在以下几方面对程序进行改进，以加

深对矩阵位移法及MATLAB语言编程的理解：

1、本程序的单元信息输入、结点坐标输入均分别隐含要求从1单元、1结点按编码的序号顺序输入，考虑如何改进使相关信息的输入可以更加灵活。

2、本程序的弹性模量仅能输入一个数值，意味着程序仅能计算由同种材料构成的结构。考虑如何改进使程序可计算由不同材料构成的组合结构。

3、本程序仅能计算两种非结点荷载类型，考虑如何编制跨中作用集中力偶、单元任意位置作用横向集中力、单元部分范围内作用均布荷载等类型的程序。

4、本程序的单元杆端力仅能计算杆端弯矩，考虑如何编制计算杆端剪力的程序

5、连续梁结构由于局部坐标系与整体坐标系相同，因此未涉及到坐标转换的问题，考虑如何编制坐标转换的程序。

6、连续梁结构由于每个结点仅有一个转角位移，所以程序中单刚、总刚、荷载矩阵等多处均采用结点号来替代结点位移号，进行下表索引和循环控制。但对于其他结构如平面刚架，每个结点一般有三个结点位移，则程序的编制过程中需要严格分清结点号和结点位移号。考虑如何编制平面刚架的计算程序。

综上所述，本章的连续梁程序体现了如何将矩阵位移法的计算方法和过程用MATLAB程序语言表达出来，重点放在程序架构和流程的建立以及算法实现方面，主要依赖手工操作－手工输入初始数据（前处理）、手工绘制计算结果（后处理），目的是使学生能够清晰、明确地掌握矩阵位移法的基本理论和概念，熟练掌握应用MATLAB语言编制、修改和调试简单程序的能力。

3.5 连续梁程序计算上机实习任务书

一、上机实习的主要内容

1、掌握连续梁计算程序输入的有关要求，应用程序进行连续梁结构的计算。

1)认真阅读连续梁计算程序的输入说明文件，进行例题的上机校验，理解并熟练

掌握程序的数据输入说明要求。

2)运用连续梁计算程序对所设计的习题进行计算。计算中，要针对每一习题的不

同要求改变相应的参数，研究参数改变对结构的影响。

3)对程序的计算结果进行分析，写出实习报告，实习报告要规范、详尽、言之有

据。

2、阅读并理解连续梁计算程序的编制思想和流程，画出程序框图，理解每一模块的功能及编程思想，然后根据自己的理解对其进行修改，使其更加完善。供参考的修改内容包括：

1)本程序的单元信息输入、结点坐标输入均分别隐含要求从1单元、1结点按编

码的序号顺序输入，考虑如何改进使相关信息的输入可以更加灵活。

2)本程序的弹性模量仅能输入一个数值，意味着程序仅能计算由同种材料构成的

结构。考虑如何改进使程序可计算由不同材料构成的组合结构。

3)本程序仅能计算两种非结点荷载类型，考虑如何编制跨中作用集中力偶、单元

任意位置作用横向集中力、单元部分范围内作用均布荷载等类型的程序。

4)本程序的单元杆端力仅能计算杆端弯矩，考虑如何编制计算杆端剪力的程序

二、上机实习的设计习题

设计多跨连续梁，通过改变下列参数，分析连续梁受力的特点：

1)改变各跨连续梁的跨度（L），探讨跨度变化对连续梁结构的受力特点；

2)改变各跨连续梁的弹性模量（E），探讨各跨采用不同材料时连续梁结构的受

力特点；

3)改变各跨连续梁的惯性矩（I），探讨各跨截面尺寸不同时连续梁结构的受力

特点；

4)改变各跨连续梁的抗弯刚度（EI），探讨各跨抗弯刚度不同时连续梁结构的

受力特点；

5)探讨跨中集中力、结点集中力偶和竖向均布荷载（如自重）的不同组合作用时

连续梁结构的受力特点；

6)改变结构的边界条件（铰支座或固定端），探讨不同边界条件时连续梁结构的

受力特点。

在以上设计中，结构的跨数、跨度、边界条件、弹性模量、惯性矩、荷载等参数均自行

设定，且应与实际工程结构相对应。当分析单一参数影响时，其他参数可取相同的设定，以利于分析单一参数对连续梁结构的影响规律。

总结以上不同要求对连续梁结构的受力影响，从而对连续梁结构的受力特点加深了解。

三、上机实习的要求

1、上机前，应在课下将每一习题的数据输入文件准备好。上机时，进行数据文件的输入、计算、修改、分析等工作，以利于节省时间。

2、上交的实习报告应包含下列内容：

1)数据输入文件及说明，包括结构单元编码、结点编码、坐标系设定、结构尺寸、

荷载等相关信息的计算简图；

2)计算结果文件及弯矩图；

3)详尽的分析；

4)总结的规律。

5)简明扼要地谈谈上机实习的体会。