2016基本平面刚架各种荷载MATLAB程序
第2章++MATLAB2016程序设计基础
2.1 MATLAB的程序界面
2.1.4当前文件夹窗口
当前文件夹(Current Folder)窗口显示了MATLAB在对文件 操作（保存、打开等）时默认的工作目录，显示的文件信 息包括文件名称、文件类型、修改日期、内容描述等。 该窗口相当于是一个资源管理器。在当前文件夹窗口中的 某一文件上单击鼠标右键，会弹出上下文菜单，可通过此 菜单实现对文件的打开、运行、重命名、复制、删除等操 作
7
2.1 MATLAB的程序界面
键
↑ ↓ ← → Ctrl+← Ctrl+ → Home
功能
重调前一行 重调下一行 左移一个字符 右移一个字符 左移一个字 右移一个字 移动到行首
8
键
End ESC Backspace PageUp PageDown Ctrl+Home Ctrl+End
功能
移动到行尾 清除一行 删除光标左侧一个字符 向前翻页 向后翻页 光标移到命令窗口首 光标移到命令窗口尾
10
2.1 MATLAB的程序界面
2.1.2 工作区窗口
工作区窗口是MATLAB的变量管理中心，每次启动 MATLAB，都会自动建立一个工作区，运行MATLAB 的程序或命令时产生的变量被加入到工作区中，除 非用特殊的命令删除某变量，否则该变量在关闭 MATLAB之前一直保存在工作区，工作区在MATLAB 运行期间一直存在，关闭MATLAB后，工作区才会自 动消除。
27
2.4 数据与数据类型
例2-3创建变量并赋值为数组。
在命令行窗口依次输入下面命令：
>>x=[1 2 3 4 5 6 7 8 9] >> y=[1,2,3;4,5,6;7,8,9] 注2: 在赋值过程中，如果变量已存在，MATLAB将使用 新值代替旧值，并以新的变量类型代替旧的变量类型。
MATLAB2016基础实例教程 第1章 MATLAB入门
1.1.4 MATLAB系统
MATLAB系统主要包括以下5个部分 桌面 数学函数库工具和开发环境 语言 图形处理 外部接口
1.2 MATLAB 2016的用户界面
MATLAB 2016的工作界面主要由标题栏、菜单栏、 工具栏、当前工作目录窗口、命令窗口、工作空 间管理窗口和历史命令窗口等组成。
1.2.5 历史窗口
选择“命令历史记录”→“停靠”命令，在显示界 面上固定显示命令历史窗口，如图所示。
在历史窗口中双击某一命令，命令窗口中将执行该 命令。
1.2.6 当前目录窗口
当前目录窗口显示如图所示，可显示或改变当前目 录，查看当前目录下的文件，单击 按钮可以在当 前目录或子目录下搜索文件。
（3）缺少步骤，未定义变量
（4）正确格式
1.3.2 功能符号
除了命令输入必须的符号外，MATLAB为了解决命 令输入过于繁琐、复杂的问题，采取了分号、续行 符及插入变量等方法。
1．分号 一般情况下，在MATLAB中命令窗口中输入命令，
则系统随机根据指令给出计算结果。若不想让 MATLAB每次都显示运算结果，只需在运算式最后加 上分号（；）。 2．续行号
1.1.3 MATLAB的特点
MATLAB的一个重要特色是它具有一系列称为工 具 箱 （ Toolbox ） 的 特 殊 应 用 子 程 序 。 工 具 箱 是 MATLAB函数的子程序库，可以分为功能性工具 箱和学科性工具箱。
所有MATLAB核心文件和各种工具箱文件都是可 读可修改的源文件，用户可通过对源程序进行修 改或加入自己编写的程序来构造新的专用工具箱。
1.2.2 功能区
MATLAB 2016将所有的功能命令分类别放置在三 个选项卡中，下面分别介绍这3个选项卡。 “主页”选项卡：单击标题栏下方的“主页” 选项卡，显示基本的“新建脚本”“新建变量” 等命令。
实验一 Matlab基本操作(2016)
实验一 MATLAB 基本操作一、实验目的1. 学习和掌握MA TLAB 的基本操作方法2. 掌握命令窗口的使用3. 熟悉MATLAB 的数据表示、基本运算二、实验内容和要求1. 实验内容1) 练习MATLAB7.0或以上版本2) 练习矩阵运算与数组运算2. 实验要求1) 每位学生独立完成，交实验报告2) 禁止玩游戏！三、实验主要软件平台装有MATLAB7.0或以上的PC 机一台四、实验方法、步骤及结果测试1. 实验方法：上机练习。
2. 实验步骤：1) 开启PC ，进入MA TLAB 。
2) 使用帮助命令，查找sqrt 函数的使用方法答: help sqrt3) 矩阵、数组运算a) 已知 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321A ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=963852741B ，求)2()(A B B A -⋅+ 答： A=[1, 2, 3; 4, 5, 6; 7, 8, 9]; B=[1, 4, 7; 2, 5, 8; 3, 6, 9]; (A+B)*(2*B-A)b) 已知⎥⎦⎤⎢⎣⎡-=33.1x ，⎥⎦⎤⎢⎣⎡=π24y ，求T xy ，y x T c) 已知⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321A ，⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=300020001B ，求A/B, A\B. d) 已知⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=987654321A ，求：(1) A 中第三列前两个元素；(2) A 中所有第二行元素；(3) A 中四个角上的元素；(4) 交换A 的第1、3列。
(5) 交换A 的第1、2行。
(6) 删除A 的第3列。
e) 已知[]321=x ，[]654=y ，求：y x *.，y x /.，y x \.，y x .^，2.^x ，x .^2。
f) 给出x=1,2,…,7时，xx sin 的值。
3）常用的数学函数a ）随机产生一个3x3的矩阵A ，求：(1) A 每一行的最大、最小值，以及最大、最小值所在的列；(2) A 每一列的最大、最小值，以及最大、最小值所在的行；(3) 整个矩阵的最大、最小值；(4) 每行元素之和；(5) 每列元素之和；(6) 每行元素之积；(7) 每列元素之积。
平面刚架电算程序说明(供阅读、熟悉程序、编程使用)
平面刚架电算程序说明一．本程序适用范围1．结构形式：刚架可具有任意几何形状，刚架结点全部为刚性结点，各杆为等截面杆。
2．支承方式：支座形式可以是固定端、铰支座、沿垂直或水平方向的棍轴支座。
3．荷载类型：荷载包括结点荷载和非结点荷载。
结点荷载包括三种：水平集中力、垂直集中力、集中力偶。
非结点荷载包括三种情况：部分均布荷载、垂直集中荷载和平行集中荷载。
二．数据符号说明1．输入基本参数：ne——单元数nj——结点数nz——支承数，为支承对应的位移分量总个数npj——结点荷载数npf——非结点荷载数e0——弹性模量nj3——位移分量总个数 nj3=nj×32．输入其它数据jm（1 to ne，1 to 2）——杆端结点码数组。
由全部NE个单元的两个端点的结点码组成。
Jm(e,1)——表示第E号单元的始端结点码Jm(e,2)——表示第E号单元的末端结点码l（1 to ne）——杆长数组an（1 to ne）——杆角数组h（1 to ne）——杆件截面高度数组b（1 to ne）——杆件截面宽度数组mj（1 to ne）——杆件截面面积数组i（1 to ne）——截面惯性距数组zc(1 to nz)——支承数组。
由全部NZ个支承对应的位移分量的编码组成。
ZC（I）——表示第I号支承对应的位移分量的编码pj(0 to npj，1 to 2)——由全部NPJ个结点荷载的数值及对应位移分量的编码组成。
当NPJ=0时，可将PJ（0，1）与PJ（0，2）设为0。
当NPJ>0时，PJ（I，1）——第I号结点荷载的数值PJ（I，2）——第I号结点荷载对应的位移分量编码pf(0 to npf，1 to 4) ——由全部NPF个非结点荷载的四种数据组成（荷载数值、位置参数、作用杆码、荷载类型码）。
当NPF=0时，可将第0行的四个元素设为0。
当NPF>0时，PF（I，1）——第I号非结点荷载的数值PF（I，2）——第I号非结点荷载的位置参数。
刚架矩阵位移法matlab计算程序
clc;%输入基本信息jd=input('请输入节点数量');free=input('请输入自由度');gj=input('请输入杆件数量');P=zeros(free,1);P=input('请输入外荷载矩阵（[x;y;z]）');d=zeros(free,1);member=zeros(1,6,gj);K=zeros(6,6,gj);k=zeros(6,6,gj)for i=1:1:gjangle(i)=input(sprintf('请输入第%d杆件角度（角度制）',i));E(i)=input(sprintf('请输入第%d杆件弹性模量（N/mm2）',i));A(i)=input(sprintf('请输入第%d杆件截面面积（mm2）',i));L(i)=input(sprintf('请输入第%d杆件长度（mm）',i));member(:,:,i)=input(sprintf('请输入第%d杆件定位向量([1 2 3 4 5 6])',i)); end%计算局部坐标系下单元刚度矩阵kfor i=1:1:gjk=(:,:,i)=E(i)*I(i)/L(i)^3*[A(i)*L(i)^2/I(i) 0 0 -A(i)*L(i)^2/I(i) 0 0;0 12 6*L(i) 0 -12 6*L(i);0 6*L(i) 4*L(i)^2 0 -6*L(i) 2*L(i)^2;-A(i)*L(i)^2/I(i) 0 0 A(i)*L(i)^2/I(i) 0 0;0 -12 -6*L(i) 0 12 -6*L(i);0 6*L(i) 2*L(i)^2 0 -6*L(i) 4*L(i)^2 ];end%计算各杆件转换矩阵TT=zeros(6,6,gj);for i=1:1:gjT=(:,:,i)=[cosd(angle(i)) sind(angle(i)) 0 0 0 0;-sind(angle(i)) cosd(angle(i)) 0 0 0 0;0 0 1 0 0 0;0 0 0 cosd(angle(i)) sind(angle(i)) 0;0 0 0 -sind(angle(i)) cosd(angle(i)) 0;0 0 0 0 0 1];end%计算整体坐标系下单元刚度矩阵KK(:,:,i)=T(:,:,i)'*k(:,:,i)*T(:,:,i);end%形成SSs=zeros(free);S=zeros(free);for i=1:1:gjfor n=1:1:6for j=1:1:6if (member(1,n,i)<=free && member(1,j,i)<=free)Ss(member(1,n,i),member(1,j,i))=K(n,j,i);endendendS=S+Ss ;Ss=zeros(free);end%计算杆间荷载作用PfFf=[FAbcosd-FSbsind;FAbsind+FSbcosd;FMb;FAecosd-FSbsind;FAesind+FSecosd;FMe]for i=1:1:gjgjl=input('请输入杆间力数量');ppp=zeros(1,4,gjl);for j=1:1:gjlppp(:,:,i)=input(sprintf('请输入第%d杆件l1,l2,W，a ([1 2 3 4])',i));l1=ppp(1,1,i);l2=ppp(1,2,i);W=ppp(1,3,i);a=ppp(1,3,i);if a==0,l1+l2==L(i)pd=1else if a==0,l1+l2<L(i)pd=2else if a>0,l1+l2==L(i)pd=3elsepf=4endendendswitch pdcase 1Qf(:,:,gjl)=[0;W*l2^2*(3*l1+l2);W*l1*l2^2;0;W*l1^2*(l1+3*l2);-W*l1^2*l2/L(i)^2]case 2Qf(:,:,gjl)=[W*cosd(a);W*sind(a)*l2^2*(3*l1+l2);W*sind(a)*l1*l2^2;W*cosd(a);W*sin(a)*l1^2*(l1+3*l2);-W*l1^2*l2/L(i)^2]case 3Qf(:,:,gjl)=[0;W*L(i)/2*(1-l1/L(i)^4*(2*L(i)^3-2*l1^2*L(i)+l1^3)-l2^3/L(i)^4*(2*L(i)-l2));W*L(i)^2/12*(1-l1^2/L(i)^4*(6*L(i)^2-8*l1*L(i)+3*l1^2)-l2^3/L(i)^4*(4*L(i)-3*l2));W*L(i)/2*(1-l1^3/L(i)^4*(2*L(i)-l1)-l2/L(i)^4*(2*L(i)^3-2*l2^2+l2^3));0;-W*L(i)^2/12*(1-l1^3/L(i)^4*(4*L(i)-3*l1)-l2^2/L(i)^4*(6*L(i)^2-8*l2*L(i)+3*l2^2))]case 4Qf(:,:,gjl)=[W*cosd(a)*(L(i)-l1-l2);W*sind(a)*L(i)/2*(1-l1/L(i)^4*(2*L(i)^3-2*l1^2*L(i)+l1^3)-l2^3/L(i)^4*(2*L(i)-l2));W*sind(a)*L(i)^2/12*(1-l1^2/L(i)^4*(6*L(i)^2-8*l1*L(i)+3*l1^2)-l2^3/L(i)^4*(4*L(i)-3*l2));W*sind(a)*L(i)/2*(1-l1^3/L(i)^4*(2*L(i)-l1)-l2/L(i)^4*(2*L(i)^3-2*l2^2+l2^3));W*sind(a)*cosd(a)*(L(i)-l1-l2);-W*sind(a)*L(i)^2/12*(1-l1^3/L(i)^4*(4*L(i)-3*l1)-l2^2/L(i)^4*(6*L(i)^2-8*l2*L(i)+3*l2^2)) ]endendFf(:,:,i)=T(:,:,i)'*Qf(:,:,i);endPf=zeros(free,1);for i=1:1:gjpfp=zeros(free,1);for j=1:1:6if member(1,j,i)<=freePfp(member(1,j,i),1)=Pf(j,1,i);endendPf=Pf+Pfp;end%P-Pf=S*dd=S\(P-Pf);%计算位移和内力v=zeros(6,1,gj);u=zeros(6,1,gj)Q=zeros(6,1,gj);F=zeros(6,1,gj);for i=1:1:gjc=1;for j=1:1:6if(member(1,j,i)<free+1)v(j,1,i)=d(c,1);c=c+1;endendendfor i=1:1:gju(:,:,i)=T(:,:,i)*v(:,:,i);endfor i=1:1:gjQ(:,:,i)=k(:,:,i)*u(:,:,i)endfor i=1:1:gjF(:,:,i)=T(:,:,i)'*Q(:,:,i);endzfl=jd*2;r=zeros(zfl,1);for i=1:1:gjfor j=1:1:4r(member(1,j,i),1)=F(j,1,i);endend%计算支座反力zfl=zfl-free;R=r(free+1:end)。
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'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
平面杆件结构计算 matlab
⎡µ1⎤ ⎡0⎤
∆1
=
⎢⎢υ1
⎥ ⎥
=
⎢⎢0⎥⎥
⎢⎣θ1 ⎥⎦ ⎢⎣0⎥⎦
⎡µ10 ⎤ ⎡0⎤
∆10
=
⎢⎢υ1
0
⎥ ⎥
=
⎢⎢0⎥⎥
⎢⎣θ10 ⎥⎦ ⎢⎣0⎥⎦
-4-
修改刚度方程有划行划列法和乘大数法，划行划列法是在原始刚度矩阵中删
去 0 位移对应的行和列，即分别划去 1 号和 10 号结点各自对应的三行和三列，
②形成局部坐标系中的单元刚度矩阵 k (e ) 和整体坐标系中的单元刚度矩阵 k (e )： 对于单元（1）：
⎡ 1.260
0
0 -1.2600
0
0⎤
⎢ ⎢
0 0.1008 0.1008
0
- 0.1008
0.1008
⎥ ⎥
k
(1
)
=
1×
109
⎢0 ⎢⎢-1.260
0.1008 0
0.1344 0
0 - 0.1008 0.0672 ⎥
i=EN(ND(id,1),1);
j=EN(ND(id,1),2);
F(((i-1)*3+1):((i-1)*3+3))=F(((i-1)*3+1):((i-1)*3+3))+EF(1:3); %等效结点力集成
F(((j-1)*3+1):((j-1)*3+3))=F(((j-1)*3+1):((j-1)*3+3))+EF(4:6); %等效结点力集成
=
(e )
kT
(e )∆(e )
−
F
(e )
E
对于单元（1）
MATLAB基本操作及环境设置
MATLAB基本操作及环境设置1.MATLAB的基本操作：-启动MATLAB：在计算机上安装MATLAB软件后，可以从开始菜单中或桌面图标启动MATLAB。
-MATLAB命令窗口：启动MATLAB后，可以看到一个命令窗口。
在命令窗口中，可以输入MATLAB命令，并执行它们。
- 基本算术操作：MATLAB可以进行基本的算术操作，如加减乘除。
例如，输入"2+3"，然后按Enter键，MATLAB将计算并显示结果。
- 变量：在MATLAB中，可以定义变量，并将值赋给它们。
例如，输入"x = 5"，然后按Enter键，MATLAB将创建变量x，并将值设为5 - 矩阵操作：MATLAB是以矩阵为基础的语言。
可以使用MATLAB的矩阵操作函数创建、修改和操作矩阵。
例如，可以使用"zeros"函数创建由0组成的矩阵，使用"eye"函数创建单位矩阵，以及使用"inv"函数计算矩阵的逆矩阵。
2.MATLAB的环境设置：- 工作目录：工作目录是MATLAB文件的位置。
可以使用"cd"命令更改工作目录。
可以使用"pwd"命令查看当前工作目录。
- 文件管理：MATLAB提供了一些函数来管理和操作文件。
可以使用"dir"函数列出当前目录中的文件和文件夹，使用"mkdir"函数创建新文件夹，使用"delete"函数删除文件等。
-图形界面：MATLAB还提供了一个图形用户界面（GUI），可以通过点击菜单和按钮来执行操作。
GUI提供了更直观和交互式的方式来使用MATLAB。
- 图形绘制：MATLAB具有强大的图形绘制功能。
可以使用"plot"函数绘制二维曲线，使用"mesh"函数绘制三维曲面等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
% 平面刚架MATLAB程序% 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2014.09 2016.03%*************************************************% 变量说明% NPOIN NELEM NVFIX NFPOIN NFPRES% 总结点数,单元数, 约束个数, 受力结点数, 非结点力数% COORD LNODS YOUNG% 结构节点坐标数组, 单元定义数组, 弹性模量% FPOIN FPRES FORCE FIXED% 结点力数组，非结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵，结点位移向量%**************************************************format short e %设定输出类型clear %清除内存变量FP1=fopen('6-6.txt','rt') %打开初始数据文件%读入控制数据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',[6,NELEM])'% 单元定义：左、右结点号，面积，惯性矩，线膨胀系数，截面高度（共计NELEM组）COORD=fscanf(FP1,'%f',[2,NPOIN])'% 坐标：x,y坐标（共计NPOIN 组）FPOIN=fscanf(FP1,'%f',[4,NFPOIN])'% 节点力（共计NFPOIN 组）：受力结点号、X方向力（向右正），% Y方向力（向上正），M力偶（逆时针正）FPRES=fscanf(FP1,'%f',[7,NFPRES])' % 均布力（共计% NFPRES 组）：单元号、荷载类型、荷载大小、距离左端长度，温差=（下端-上端）梯形上边。
下边(改)% 荷载类型1-均布荷载2-横向集中力3-纵向集中力4-三角形荷载5-温度荷载6-梯形荷载FIXED=fscan f(FP1,'%f',NVFIX)'% 约束信息：约束对应的位移编码（共计NVFIX 组）%---------------------------------------------------------HK=zeros(3*NPOIN,3*NPOIN); % 张成总刚矩阵并清零FORCE=zeros(3*NPOIN,1); % 张成总荷载向量并清零%形成总刚for i=1:NELEM % 对单元个数循环% 生成局部单刚(局部坐标) 右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD);% 坐标转换矩阵EKT=T’*EK*T；% 生成整体单刚(整体坐标系)% 组成总刚按3*3子块加入总刚中（共计4块)for j=1:2 %对行进行循环---按结点号循环N1=LNODS(i,j)*3; % j结点第3个位移的整体编码for k=1:2 %对列进行循环---按结点号循环N2=LNODS(i,k)*3; % k结点第3个位移的整体编码HK((N1-2):N1,(N2-2):N2)=HK((N1-2):N1,(N2-2):N2)...+EKT(j*3-2:j*3,k*3-2:k*3);% 单刚3 x 3子块叠加到总刚中endend% 由结点力与非结点力生成总荷载向量列阵for i=1:NFPOIN % 对结点荷载个数进行循环N1=FPOIN(i,1); % 作用荷载的结点号N1=N1*3-3; % 该结点号对应第一个位移编码- 1for j=1:3FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取结点荷载endend% 计算由非结点荷载引起的等效结点荷载for i=1:NFPRES % 对非结点荷载个数进行循环F0=ele_FPRES(i,FPRES,LNODS,COORD,YOUNG); %计算单元固端力% 对单元局部杆端力要进行坐标转换ele=FPRES(i,1); % 取荷载所在的单元号T=zbzh(ele,LNODS,COORD); % 坐标转换矩阵F0=T’*F0;NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号% 将单元固端力变成等效结点荷载（注意固端力与等效结点荷载符号相反）FORCE((3*NL-2):3*NL)=FORCE((3*NL-2):3*NL)-F0(1:3);FORCE((3*NR-2):3*NR)=FORCE((3*NR-2):3*NR)-F0(4:6);end% 总刚、总荷载进行边界条件处理for j=1:NVFIX % 对约束个数进行循环N1=FIXED(j);HK(1:3*NPOIN,N1)=0; HK(N1,1:3*NPOIN)=0; HK(N1,N1)=1;% 将零位移约束对应的行、列变成零，主元变成1FORCE(N1)=0;end%---------------------------------------------------------DISP=HK\FORCE % 方程求解，HK先求逆再与力向量左乘% 求结构各个单元内力EDISP=zeros(6,1); % 单元位移列向量清零for i=1:NELEM % 对单元个数进行循环for j=1:2 %对杆端循环% i单元左右端结点号*3 = 该结点的最后一个位移编码N1=LNODS(i,j)*3;% 取一端的单元位移列向量EDISP(3*j-2:3*j)=DISP(N1-2:N1);end% 生成局部单刚(局部坐标) 右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD); % 坐标转换矩阵FE=EK*T*EDISP; %计算局部坐标杆端力（由结点位移产生）for j=1:NFPRESif FPRES(j,1) == i %成立时，当前单元上有非结点荷载F0=ele_FPRES(j,FPRES,LNODS,COORD,YOUNG);%单元固端力FE=FE+F0; % 考虑由非结点荷载引起的杆端力endendFE % 打印杆端力endend%-------------------------------------------------------ele_FPRES.m %计算单元固端力函数（正方向：X向右Y向上M逆时针）% 入口参数：荷载序号，荷载信息，单元信息，结点坐标% 出口参数：单元固端力——左右两端的轴力、剪力、弯矩function F0=ele_FPRES(iFPRES,FPRES,LNODS,COORD, E)ele=FPRES(iFPRES,1); %取荷载所在的单元号G=FPRES(iFPRES,3); %单元荷载大小C=FPRES(iFPRES,4); %单元荷载与左端距离W= FPRES(iFPRES,5); %单元下上端温差S= FPRES(iFPRES,6); %梯形长边荷载X= FPRES(iFPRES,6); %梯形短边荷载H= LNODS(ele,6);%单元截面高度P = LNODS(ele,5);%单元截面线膨胀系数NL=LNODS(ele,1); NR=LNODS(ele,2); %单元的左右结点号dx=COORD(NR,1)-COORD(NL,1); % x 坐标差dy=COORD(NR,2)-COORD(NL,2); % y 坐标差L=sqrt(dx^2+dy^2); %单元长度% 计算公式中一些常出现的项D=L-C; C1=C/L; C2=C1*C1; C3=C1*C2;B1=D/L; B2=B1/L;F0=[0;0;0;0;0;0]; %单元固端力清零switch FPRES(iFPRES,2)case 1 %均布荷载F0(2)=-G*C*(2-2*C2+C3)/2.0;F0(3)=-G*C*C*(6-8*C1+3*C2)/12.0;F0(5)=-G*C-F0(2);F0(6)=G*C*C*C1*(4-3*C1)/12.0;case 2 %横向集中力F0(2)=-G*B1*B2*(L+2*C);F0(3)=-G*C*B1*B1;F0(5)=-G*C2*(L+2*D)/L;F0(6)=G*D*C2;case 3 %纵向集中力F0(1)=-G*B1;F0(4)=-G*C1;case 4 %三角形荷载F0(2)=-7*G*L/20;F0(3)=G*L^2/20;F0(5)=3*G*L/20;F0(6)=G*L*L/30;case 5 %温度荷载F0(3)=E*I*W*P/H;F0(6)=-E*I*W*P/H;case 6 %梯形荷载（改）F1(2)=-X*C*(2-2*C2+C3)/2.0;F1(3)=-X*C*C*(6-8*C1+3*C2)/12.0;F2(5)=-X*C-F0(2);F2(6)=X*C*C*C1*(4-3*C1)/12.0;F3(2)=-7*(S-G)*L/20;F3(3)=(S-G)*L^2/20;F4(5)=3*(S-G)*L/20;F4(6)=(S-G)*L*L/30;F0(2)= F1(2)+ F3(2);F0(3)= F1(3)+ F3(3);F0(5)= F2(5)+ F4(5);F0(6)= F2(6)+ F4(6);endreturnele_EK.m % 计算单元刚度矩阵函数EK% 入口参数：单元号、单元信息数组、结点坐标、弹性模量% 出口参数：局部单元刚度矩阵EKfunction EK=ele_EK(i,LNODS,COORD,E)NL=LNODS(i,1); NR=LNODS(i,2); %左右结点号dx=COORD(NR,1)-COORD(NL,1); % x 坐标差dy=COORD(NR,2)-COORD(NL,2); % y 坐标差L=sqrt(dx^2+dy^2); %单元长度A=LNODS(i,3); I=LNODS(i,4); %面积；惯性矩% 生成单刚(局部坐标) 右手坐标系EK =[E*A/L 0 0 -E*A/L 0 0;...0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2;...0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;...-E*A/L 0 0 E*A/L 0 0;...0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;...0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L]; return%--------------------------------------------------------zbzh.m % 形成第i单元的坐标转换矩阵函数T(6,6)% 入口参数：单元号，单元信息，结点坐标% 出口参数：坐标转换矩阵（整体向局部投影）function T=zbzh(i,LNODS,COORD)NL=LNODS(i,1); %左结点号NR=LNODS(i,2); %右结点号dx=COORD(NR,1)-COORD(NL,1); % x 坐标差dy=COORD(NR,2)-COORD(NL,2); % y 坐标差L=sqrt(dx^2+dy^2); % 单元长度c=dx/L; % cos a (与x 轴夹角余弦）s=dy/L; % sin aT=[ c s 0 0 0 0;...-s c 0 0 0 0;...0 0 1 0 0 0;...0 0 0 c s 0;...0 0 0 -s c 0;...0 0 0 0 0 1];return。