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

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

入结构力学领域中而产生的一种方法，而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,NVF IX)

●说明：

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 列的主对角线元素HK(N1,N1)乘以一个大数，如10

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)