下载文档原格式
实验一 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号非结点荷载的位置参数。
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语言编制连续梁结构的计算程序。
矩阵位移法的解题思路是将结构离散为单元（杆件），建立单元杆端力与杆端位移之间的关系－单元刚度方程；再将各单元集成为原结构，在满足变形连续条件和平衡条件时，建立整体刚度方程；在边界条件处理完毕后，由整体刚度方程解出节点位移，进而求出结构内力。
用矩阵位移法计算连续梁的步骤如下：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.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"函数绘制三维曲面等。
matlab有限元编程荷载-概述说明以及解释1.引言1.1 概述概述部分的内容可以包括以下内容：有限元方法是一种广泛应用于工程领域的数值计算方法，它通过将复杂的连续体结构分割为一系列简单的子结构，然后利用数学方法为每个子结构建立相应的数学模型，从而得到整个结构的行为特性。
这种分割的过程通常被称为网格划分，而每个子结构则称为有限元。
MATLAB作为一种功能强大的数学软件工具，被广泛应用于有限元方法的编程与分析。
MATLAB提供了大量的工具箱和函数，可以方便地实现有限元方法的各个步骤，包括网格划分、单元构造、边界条件的施加以及结果的可视化分析等。
本文将重点介绍MATLAB在有限元编程中的应用，特别是在荷载分析方面的应用。
荷载分析是有限元分析的核心内容之一，它通过施加不同的荷载条件，分析结构在荷载作用下的变形、位移和应力等。
荷载分析的准确性对于工程设计以及结构安全性的评估至关重要。
文章将首先介绍有限元方法的基本原理，包括对结构的离散化、单元的建立和组装，以及求解过程中的矩阵运算等。
然后，针对荷载分析的特点，将详细介绍MATLAB有限元编程中的荷载处理方法，包括荷载的施加、荷载类型的选择以及荷载与结构响应的耦合关系等。
通过本文的学习，读者可以了解到MATLAB在有限元编程中的应用，并且掌握荷载分析的基本原理和方法。
同时，也可以对MATLAB有限元编程在实际工程问题中的应用进行进一步的探索和研究。
总而言之，本文将为读者提供一个全面而系统的MATLAB有限元编程荷载分析的引导，帮助读者理解有限元方法的基本原理和应用，提高工程设计和结构安全性评估的能力。
1.2 文章结构本文旨在介绍MATLAB在有限元编程中的应用，重点讲解荷载分析的基本原理以及MATLAB有限元编程中的荷载处理方法。
该文章分为引言、正文和结论三个部分。
引言部分对文章进行了概述，包括本文的目的和总结。
在概述中，我们会介绍有限元方法的简要背景和意义，以及MATLAB在该领域中的重要性。
实验项目及学时安排实验一 MATLAB环境的熟悉与基本运算 2学时实验二 MATLAB数值计算实验 2学时实验三 MATLAB数组应用实验 2学时实验四 MATLAB符号计算实验 2学时实验五 MATLAB的图形绘制实验 2学时实验六 MATLAB的程序设计实验 2学时实验七 MATLAB工具箱Simulink的应用实验 2学时实验八 MATLAB图形用户接口GUI的应用实验 2学时实验一 MATLAB环境的熟悉与基本运算一、实验目的1．熟悉MATLAB开发环境2．掌握矩阵、变量、表达式的各种基本运算二、实验基本知识1.熟悉MATLAB环境:MATLAB桌面和命令窗口、命令历史窗口、帮助信息浏览器、工作空间浏览器、文件和搜索路径浏览器。
2.掌握MATLAB常用命令3.MATLAB变量与运算符变量命名规则如下：（1）变量名可以由英语字母、数字和下划线组成（2）变量名应以英文字母开头（3）长度不大于31个（4）区分大小写MATLAB中设置了一些特殊的变量与常量，列于下表。
MATLAB运算符，通过下面几个表来说明MATLAB的各种常用运算符4.MATLAB的一维、二维数组的寻访表6 子数组访问与赋值常用的相关指令格式5.MATLAB的基本运算表7 两种运算指令形式和实质涵的异同表6.MATLAB的常用函数表8 标准数组生成函数表9 数组操作函数三、实验容1、学习使用help命令，例如在命令窗口输入help eye，然后根据帮助说明，学习使用指令eye（其它不会用的指令，依照此方法类推）2、学习使用clc、clear，观察command window、command history和workspace等窗口的变化结果。
3、初步程序的编写练习，新建M-file，保存（自己设定文件名，例如exerc1、exerc2、 exerc3……），学习使用MATLAB的基本运算符、数组寻访指令、标准数组生成函数和数组操作函数。
% 平面刚架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。
