平面三角形单元Matlab程序

%变量说明

%NPOIN NELEM NVFIX NFORCE NNODE

%总结点数，单元数，受约束边界点数，结点力数, 单元结点数

%COORD LNODS YOUNG POISS THICK %结构结点坐标数组，单元定义数组，弹性模量，泊松比，厚度

%FORCE FIXED BMA TX DMATX SMATX %节点力数组，约束信息数组，单元应变矩阵，单元弹性矩阵，单元应力矩阵%AREA HK ASLOD ASDISP FP1

%单元面积，总体刚度矩阵，总体荷载向量，结点位移向量，数据文件指针format short e

clear

FP1=fopen('C:\Users\Administrator\Desktop\input.txt','rt');

NPION=fscanf(FP1,'%d',1); %结点个数（结点编码总数）NELEM=fscanf(FP1,'%d',1); %单元个数（单元编码总数）NFORCE=fscanf(FP1,'%d',1); %结点荷载个数

NVFIX=fscanf(FP1,'%d',1); %受约束边界点数

YOUNG=fscanf(FP1,'%e',1); %弹性模量

POISS=fscanf(FP1,'%f',1); %泊松比

THICK=fscanf(FP1,'%d',1); %厚度

LNODS=fscanf(FP1,'%d',[3,NELEM])'; %单元定义数组（单元结点号）COORD=fscanf(FP1,'%f',[2,NPION])'; %结点坐标数组

FORCE=fscanf(FP1,'%f',[3,NFORCE])'; %结点力数组

FIXED=fscanf(FP1,'%d',[3,inf])'; %约束信息数组

%引用所需的全局变量

%global NPION NELEM COORD LNODS YOUNG POISS

%global BMA TX DMATX SMA TX AREA

%生成弹性矩阵D

a=YOUNG/(1-POISS^2);

DMATX(1,1)=1*a;

DMATX(1,2)=POISS*a;

DMATX(2,1)=POISS*a;

DMATX(2,2)=1*a;

DMATX(3,3)=(1-POISS)*a/2;

for i=1:NELEM; %i为当前所计算的单元号

%计算当前单元的面积

AREA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);...

1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);...

1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);])/2;

end

%生成应变矩阵B

for j=0:2

b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);

c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);

end

BMA TX=[b(1) 0 b(2) 0 b(3) 0;...

0 c(1) 0 c(2) 0 c(3);...

c(1) b(1) c(2) b(2) c(3) b(3)]/(2*AREA);

SMA TX=DMATX*BMATX; %求应力矩阵S=D*B

% 所引用的全局变量：

%global NPION NELEM NVFIX LNODS ASTIF THICK

%global BMA TX SMA TX AREA FIXED

HK=seros(2*NPION,2*NPION); %张成特定大小总刚矩阵并置0

for i=1:NELEM

EK=BMA TX'*SMATX*THICK*AREA; %求解单元刚度矩阵

a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号

for j=1:3

for k=1:3

HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+EK(j*2-1:j*2,k* 2-1:k*2);

%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中

end

end

end

%将约束信息加入总刚（置0置1法）

NUM=1; %计数器(当前已分析的节点数)

i=0; %计数器(当前已处理的约束数)

tmp(NVFIX)=0; %临时存被约束的序号

while i

for j=-1:0

if FIXED(NUM,j+3)==1 %若发现约束

i=i+1; %计数器自增

tmp(i)=FIXED(NUM)*2+j; %求约束序号

end

end

NUM=NUM+1;

end

for i=1:NVFIX

HK(1:2*NPION,tmp(i))=0; %将一约束序号处的总刚列向量清0

HK(tmp(i),1:2*NPION)=0; %.将一约束序号处的总刚行向量清0

HK(tmp(i),tmp(i))=1; %将行列交叉处的元素置为1

end

%所需引用的全局变量

%global ASLOD NPION NFORCE FORCE

ASLOD(1:2*NPION)=0; %张成特定大小的0向量

for i=1:NFORCE

ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);

end

%将有约束处的荷载置为0

NUM=1;

i=0;

tmp(NVFIX)=0;

while i

for j=-1:0

if FIXED(NUM,j+3)==1

i=i+1;

tmp(i)=FIXED(NUM)*2+j;

end

end

NUM=NUM+1;

end

for i=1:NVFIX

ASLOD(tmp(i))=0;

end

ASDISP=HK\ASLOD'; %计算结点位移向量

%所引用的全局变量

%global NELEM NPION SMA TX ASDISP LNODS

ELEDISP(1:6)=0; %当前单元的结点位移向量for i=1:NELEM

for j=1:3

ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);

end

STRESS=SMA TX*ELEDISP'; %求内力end

平面三角形单元有限元程序设计

. 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m ，E=200GPa ，=0.25，t=0.1m ，忽略自重。试计算薄板的位移及应力分布。 要求： 1. 编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2. 采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3. 提交程序编写过程的详细报告及计算机程序； 4. 所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则

刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算 单元刚度的集成：[ ][][][][][]' '''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +?++=? =?==?==?=?????? 边界约束的处理：划0置1法 X Y P X Y P

matlab有限元分析实例

MATLAB: MATLAB是美国MathWorks公司出品的商业数学软件，用于数据分析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人，控制系统等领域。 MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室），软件主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。 MATLAB有限元分析与应用:

《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书，作者是卡坦，译者是韩来彬。 内容简介： 《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用，书中的每个示例都以交互的方式求解，使读者很容易就能把MATLAB用于有限分析和应用。另外，《MATLAB有限元分析与应用》还提供了大量免费资源。 《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。《MATLAB有限元分析与应用》由简单到复杂，循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题，具有很高的工程应用价值。

平面三角形单元常应变单元matlab程序的编制

三角形常应变单元程序的编制与使用 有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。 有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角 形有更高的计算精度。 Matlab语言是进行矩阵运算的强大工具，因 此，用Matlab语言编写有限元中平面问题的程序 有优越性。本章将详细介绍如何利用Matlab语言 编制三角形常应变单元的计算程序，程序流程图见 图1。 有限元法中三节点三角形分析结构的步骤如 下： 1）整理原始数据，如材料性质、荷载条件、约 束条件等，离散结构并进行单元编码、结点 编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩 阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。计算结果整理包括位移和应 力两个方面；位移计算结果一般不需要特别 的处理，利用计算出的节点位移分量，就可 画出结构任意方向的位移云图；而应力解的 误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。图1 程序流程图

1.1 程序说明 %******************************************************************* % 三角形常应变单元求解结构主程序 %******************************************************************* ●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。 ●基本思想：单元结点按右手法则顺序编号。 ●荷载类型：可计算结点荷载。 ●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算 法的全过程，即实现程序流程图的程序表达。 %----------------------------------------------------------------------------------------------------- 1 程序准备 format short e %设定输出类型 clear all %清除所有已定义变量 clc %清屏 ●说明： format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法； clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况； clc －清屏，使屏幕在本程序运行开始时 %----------------------------------------------------------------------------------------------------- 2 全局变量定义 global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICK global FORCE FIXED global BMATX DMATX SMATX AREA global ASTIF ASLOD ASDISP global FP1 ●说明： NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度

主成分分析报告matlab程序

Matlab编程实现主成分分析 .程序结构及函数作用 在软件Matlab中实现主成分分析可以采取两种方式实现：一是通过编程来实现；二是直接调用Matlab种自带程序实现。下面主要主要介绍利用Matlab的矩阵计算功能编程实现主成分分析。 1程序结构 2函数作用 Cwstd.m——用总和标准化法标准化矩阵 Cwfac.m——计算相关系数矩阵；计算特征值和特征向量；对主成分进行排序；计算各特征值贡献率；挑选主成分（累计贡献率大于85%），输出主成分个数；计算主成分载荷 Cwscore.m——计算各主成分得分、综合得分并排序 Cwprint.m——读入数据文件；调用以上三个函数并输出结果

3.源程序 3.1 cwstd.m总和标准化法标准化矩阵 %cwstd.m,用总和标准化法标准化矩阵 function std=cwstd(vector) cwsum=sum(vector,1); %对列求和 [a,b]=size(vector); %矩阵大小,a为行数,b为列数 for i=1:a for j=1:b std(i,j)= vector(i,j)/cwsum(j); end end 3.2 cwfac.m计算相关系数矩阵 %cwfac.m function result=cwfac(vector); fprintf('相关系数矩阵:\n') std=CORRCOEF(vector) %计算相关系数矩阵 fprintf('特征向量(vec)及特征值(val)：\n') [vec,val]=eig(std) %求特征值(val)及特征向量(vec) newval=diag(val) ; [y,i]=sort(newval) ; %对特征根进行排序，y为排序结果，i为索引fprintf('特征根排序：\n') for z=1:length(y) newy(z)=y(length(y)+1-z); end fprintf('%g\n',newy) rate=y/sum(y); fprintf('\n贡献率：\n') newrate=newy/sum(newy) sumrate=0; newi=[]; for k=length(y):-1:1 sumrate=sumrate+rate(k); newi(length(y)+1-k)=i(k); if sumrate>0.85 break; end end %记下累积贡献率大85%的特征值的序号放入newi中fprintf('主成分数：%g\n\n',length(newi)); fprintf('主成分载荷：\n') for p=1:length(newi)

层次分析报告法及matlab程序

层次分析法建模 层次分析法（AHP－Analytic Hierachy process）---- 多目标决策方法 70 年代由美国运筹学家T·L·Satty提出的，是一种定性与定量分析相结合的多目标决策分析方法论。吸收利用行为科学的特点，是将决策者的经验判断给予量化，对目标（因素）结构复杂而且缺乏必要的数据情况下，採用此方法较为实用，是一种系统科学中，常用的一种系统分析方法，因而成为系统分析的数学工具之一。 传统的常用的研究自然科学和社会科学的方法有： 机理分析方法：利用经典的数学工具分析观察的因果关系； 统计分析方法：利用大量观测数据寻求统计规律，用随机数学方法描述（自然现象、 社会现象）现象的规律。 基本内容：（1）多目标决策问题举例AHP建模方法 （2）AHP建模方法基本步骤 （3）AHP建模方法基本算法 （3）AHP建模方法理论算法应用的若干问题。 参考书：1、姜启源，数学模型（第二版，第9章；第三版，第8章），高等教育出版社 2、程理民等，运筹学模型与方法教程，（第10章），清华大学出版社 3、《运筹学》编写组，运筹学（修订版），第11章，第7节，清华大学出版社 一、问题举例： A．大学毕业生就业选择问题 获得大学毕业学位的毕业生，“双向选择”时，用人单位与毕业生都有各自的选择标准和要求。就毕业生来说选择单位的标准和要求是多方面的，例如： ①能发挥自己的才干为国家作出较好贡献（即工作岗位适合发挥专长）； ②工作收入较好（待遇好）； ③生活环境好（大城市、气候等工作条件等）； ④单位名声好（声誉-Reputation）； ⑤工作环境好（人际关系和谐等） ⑥发展晋升（promote, promotion）机会多（如新单位或单位发展有后劲）等。 问题：现在有多个用人单位可供他选择，因此，他面临多种选择和决策，问题是他将如何作出决策和选择？——或者说他将用什么方法将可供选择的工作单位排序？

平面三角形单元

第八章 平面问题的有限元分析及三角形单元的应用 第一节 概述 分析弹性力学平面问题时，最简单的单元式由三个结点组成的三角形单元。当用以分析平面应力问题时，可将其视为三角板；当用以分析平面应变问题时，则可式为三棱柱。各单元在结点处为铰结。图8-1所示位移悬臂梁离散为三角形单元的组合体 以矩阵形式列出弹性力学平面问题的基本量和基本方程。 谈形体所受体力分量可表示为 [ ] T y x y x p p p p p =??? ? ????= (8-1) 所受面力分量可表示为 [ ] T y x y x p p p p p =??? ? ????= （8-2） 体内任一点应力分量可表示为 []T xy y x τδδδ= （8-3） 任一点的应变分量可表示为 []T xy y x γεεε= （8-4） 任一点的位移分量可表示为 []T v u =δ （8-5） 弹性力学平面问题的几何方程的矩阵表达式为 ?? ???? ???????? ??? ???+??????=????????????=x u y v y v x u xy y x εεεε （8-6） 平面应力问题的物理方程的矩阵表达式为 ? ???? ? ?????????? ????????? ?--= ????? ? ??????xy y x xy y x E γεεμμμ μτσσ210 0010112 （8-7） 或简写成为 εσD = （8-8） 式中

???? ?? ? ?????? ?--=210 0010112μμμ μ E D （8-9） 称为弹性矩阵。 平面应变问题的物理方程也可写成式（8-8），但须将式（8-9）中的E 换成 2 1μ -E ，μ换成 2 1μμ -，因此得出 ???? ?? ????????? ?? ?-----+-= )1(2210 00110 11)21)(1()1(2 2 μμμμμμ μμμE D （8-10） 平衡微分方程及边界条件也可以用矩阵表示，但弹性力学有限元位移法中，通常用虚功 方程代替平衡微分方程和应力边界条件。虚功方程的矩阵表达式为 ?????***=+tdxdy tds p f ptdxdy f T T σε （8-11） 式中：[ ] T v u f ** * =，表示虚位移； []T xy x x * ***=γεεε，表示与虚位移相对应的虚应变。 为了便于计算，作用于弹性体上的体力和面力替换为作用在结点上的集中力，即等效结 点荷载。设作用于各个结点上的外力分量用如下列阵来表示 []T n n V U V U V U F ?=2211 与这些结点外力分量相对应得结点虚位移分量列阵为 []T n n v u v u v u * ******?=2211δ 则外力在虚位移上做的虚功为 F v V u U v V u U v V u U T n n n n ** *****=++?++++δ22221111 如平面弹性体的厚度为t ，该虚功除以t ，即可得出单位厚度薄板上的外力虚功。于是，式（8-11）所示虚功方程可写成 ??**=tdxdy F T T σεδ （8-11） 虚功方程不仅仅应用于弹性力学，也可用于塑性力学。其应用条件是：只要变形体的全部外力和应力满足平衡方程；位移是微小的，并满足边界条件，位移与应变满足几何方程。

有限元的MATLAB解法

有限元的MATLAB解法 1.打开MATLAB。 2.输入“pdetool”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是partial differential equations的缩写），需要的话可点击Options菜单下Grid命令，打开栅格。 3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型，可画矩形R，椭圆E，圆C，然后在Set formula栏进行编辑并（如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标） 用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。 4.用左键双击矩形进行坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。 5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点击

“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。 6.进入PDE模式：点击"PDE"菜单下“PDE Mode”命令，进入PDE 模式，单击“PDE Specification”,设置方程类型，“Elliptic”为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。 7.对模型进行剖分：点击“Mesh”中“Initialize Mesh”进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。 8.进行计算：点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解，u值即为Hz或者Ez。 9.单击“Plot”菜单下“Parameters”选项，打开“Plot Selection”对话框。选中Color,Height(3-D plot)和Show mesh三项，然后单击“Plot”按钮，显示三维图形解。 10.如果要画等值线图和矢量场图，单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Contour和Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。 11.将计算结果条件和边界导入MATLAB中：点击“Export Solution”,再点击“Mesh”中“Export Mesh”。

平面三角形单元Matlab程序

%变量说明 %NPOIN NELEM NVFIX NFORCE NNODE %总结点数，单元数，受约束边界点数，结点力数, 单元结点数 %COORD LNODS YOUNG POISS THICK %结构结点坐标数组，单元定义数组，弹性模量，泊松比，厚度 %FORCE FIXED BMA TX DMATX SMATX %节点力数组，约束信息数组，单元应变矩阵，单元弹性矩阵，单元应力矩阵%AREA HK ASLOD ASDISP FP1 %单元面积，总体刚度矩阵，总体荷载向量，结点位移向量，数据文件指针format short e clear FP1=fopen('C:\Users\Administrator\Desktop\input.txt','rt'); NPION=fscanf(FP1,'%d',1); %结点个数（结点编码总数）NELEM=fscanf(FP1,'%d',1); %单元个数（单元编码总数）NFORCE=fscanf(FP1,'%d',1); %结点荷载个数 NVFIX=fscanf(FP1,'%d',1); %受约束边界点数 YOUNG=fscanf(FP1,'%e',1); %弹性模量 POISS=fscanf(FP1,'%f',1); %泊松比 THICK=fscanf(FP1,'%d',1); %厚度 LNODS=fscanf(FP1,'%d',[3,NELEM])'; %单元定义数组（单元结点号）COORD=fscanf(FP1,'%f',[2,NPION])'; %结点坐标数组 FORCE=fscanf(FP1,'%f',[3,NFORCE])'; %结点力数组 FIXED=fscanf(FP1,'%d',[3,inf])'; %约束信息数组 %引用所需的全局变量 %global NPION NELEM COORD LNODS YOUNG POISS %global BMA TX DMATX SMA TX AREA %生成弹性矩阵D a=YOUNG/(1-POISS^2); DMATX(1,1)=1*a; DMATX(1,2)=POISS*a; DMATX(2,1)=POISS*a; DMATX(2,2)=1*a; DMATX(3,3)=(1-POISS)*a/2; for i=1:NELEM; %i为当前所计算的单元号 %计算当前单元的面积 AREA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);... 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);... 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);])/2; end %生成应变矩阵B

基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板

基于Matlab语言的按平面三角形单元划分的结构有限元程序设计 专业：建筑与土木工程 班级：建工研12-2 姓名：韩志强 学号： 471220580

基于Matlab语言的按平面三角形单元划分 结构有限元程序设计 一、有限单元发及Matlab语言概述 1. 有限单元法 随着现代工业、生产技术的发展，不断要求设计高质量、高水平的大型、复杂和精密的机械及工程结构。为此目的，人们必须预先通过有效的计算手段，确切的预测即将诞生的机械和工程结构，在未来工作时所发生的应力、应变和位移因此，需要寻求一种简单而又精确的数值分析方法。有限单元法正是适应这种要求而产生和发展起来的一种十分有效的数值计算方法。 有限元法把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。 有限单元法基本步骤如下： (1)结构离散：结构离散就是建立结构的有限元模型，又称为网格划分或单元划分，即将结构离散为由有限个单元组成的有限元模型。在该步骤中，需要根据结构的几何特性、载荷情况等确定单元体内任意一点的位移插值函数。 (2)单元分析：根据弹性力学的几何方程以及物理方程确定单元的刚度矩阵。 (3)整体分析:把各个单元按原来的结构重新连接起来，并在单元刚度矩阵的基础上确定结构的总刚度矩阵，形成如下式所示的整体有限元线性方程： {}[]{}δ F=① K 式中，{}F是载荷矩阵,[]K是整体结构的刚度矩阵,{}δ是节点位移矩阵。 (4)载荷移置:根据静力等效原理，将载荷移置到相应的节点上，形成节点载荷矩阵。 (5)边界条件处理:对式①所示的有限元线性方程进行边界条件处理。 (6)求解线性方程:求解式①所示的有限元线性方程，得到节点的位移。在该步骤中，若有限元模型的节点越多，则线性方程的数量就越多，随之有限元分析的计算量也将越大。 (7)求解单元应力及应变根据求出的节点位移求解单元的应力和应变。

matlab动力学分析程序详解

1 1.微分方程的定义 对于duffing 方程03 2 =++x x x ω ，先将方程写作??? --==3 1122 21x x x x x ω function dy=duffing(t,x) omega=1;%定义参数 f1=x(2); f2=-omega^2*x(1)-x(1)^3; dy=[f1;f2]; 2.微分方程的求解 function solve (tstop) tstop=500;%定义时间长度 y0=[0.01;0];%定义初始条件 [t,y]=ode45('duffing',tstop,y0,[]); function solve (tstop) step=0.01;%定义步长 y0=rand(1,2);%随机初始条件 tspan=[0:step:500];%定义时间范围 [t,y]=ode45('duffing',tspan,y0); 3.时间历程的绘制 时间历程横轴为t ，纵轴为y ，绘制时只取稳态部分。 plot(t,y(:,1));%绘制y 的时间历程 xlabel('t')%横轴为t ylabel('y')%纵轴为y grid;%显示网格线

2 axis([460 500 -Inf Inf])%图形显示范围设置 4.相图的绘制 相图的横轴为y ，纵轴为dy/dt ，绘制时也只取稳态部分。红色部分表示只取最后1000个点。 plot(y(end-1000:end ,1),y(end-1000:end ,2));%绘制y 的时间历程 xlabel('y')%横轴为y ylabel('dy/dt')%纵轴为dy/dt grid;%显示网格线 5.Poincare 映射的绘制 对于不同的系统，Poincare 截面的选取方法也不同 对于自治系统一般每过其对应线性系统的固有周期，截取一次 对于非自治系统，一般每过其激励的周期，截取一次 例程：duffing 方程03 2=++x x x ω 的poincare 映射 function poincare(tstop) global omega; omega=1; T=2*pi/omega;%线性系统的周期或激励的周期 step=T/100;%定义步长为T/100 y0=[0.01;0];%初始条件 tspan=[0:step:100*T];%定义时间范围 [t,y]=ode45('duffing',tspan,y0); for i=5000:100:10000%稳态过程每个周期取一个点 plot(y(i,1),y(i,2),'b.'); hold on;% 保留上一次的图形 end xlabel('y');ylabel('dy/dt');

平面三角形单元有限元程序设计

平面三角形单元有限元程序设计 P 9 m 9 m 一、题目 如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。已知:P=150N/m,E=200GPa,=0、25,t=0、1m,忽略自重。试计算薄板的位移及应力分布。 要求: 1.编写有限元计算机程序,计算节 点位移及单元应力。(划分三角形单元,单元数不得少于30个); 2.采用有限元软件分析该问题(有 限元软件网格与程序设计网格必须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值); 3.提交程序编写过程的详细报告及计算机程序; 4.所有同学参加答辩,并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下: 1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2)单元分析,建立单元刚度矩阵。 3)整体分析,建立总刚矩阵。 4)建立整体结构的等效节点荷载与总荷载矩阵 5)边界条件处理。 6)解方程,求出节点位移。 7)求出各单元的单元应力。 8)计算结果整理。

一、程序设计 网格划分 如图,将薄板如图划分为6行,并建立坐标系,则 刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算:(1)每个单元对应2种单元刚度矩阵中的哪一种; (2)该单元对应总刚度矩阵的那几行哪几列 (3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环:(1)最外层:逐行计算 (2)中间层:该行逐个计算 (3)最里层:区分为第 奇/偶 数个计算 X Y P X Y P 节点编号 单元编号

平面三角形单元有限元程序的设计说明

. . P 9 m 9 m 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。试计算薄板的位移及应力分布。 要求： 1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3.提交程序编写过程的详细报告及计算机程序； 4.所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则

刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算 单元刚度的集成：[ ][][][][][]' '''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +?++=? =?==?==?=?????? 边界约束的处理：划0置1法 X Y P X Y P 节点编号 单元编号

平面三角形与空间四面体之间的类比

平面三角形与空间四面体之间的类比 “类比是伟大的引路人，求解立体几何问题往往有赖于平面几何中的类比问题”（波利亚）。新教材中引入类比这一内容，从根本上改变了我以往对数学的看法。虽然我以前也知道到类比，但却不敢把它作为一种数学方法理直气壮地在课堂上讲授，让学生使用。如今总算可以放开手脚，大胆应用了。 首先，平面三角形是平面几何中的一个基本图形，而四面体是立体几何中的一个基本图形。二者之间有着密切的联系，同时它们之间的联系体现了平面与空间的联系，一维空间与二维空间的联系，进一步可能有助于对多维空间的理解。 一、从概念上看：三角形是边数最少的多边形，四面体是面数最少的多面体。 二、三角形的任意两边之和大于第三边。四面体任意三个面的面积之和大于第四个面的面积。 三、任意一个三角形都有一个外接圆，即不共线三点确定一个圆，这个圆圆心称为三角形的外心，外心是各边垂直平分线的交点，外心到三角形各顶点距离相等。任意一个四面体都有一个外接球，即不共面四点确定一个球；这个球的球心在四面体各个面内的射影是各个面的外心，且它到四面体各顶点的距离也相等。 四、任意一个三角形都有一个内切圆，圆心称为三角形的内心，内心到各边距离相等，是三内角平分线的交点； 且设三角形的周长为c,内切圆半径为r,则三角形的面积为。任意一个四面体都有一个内切球，球心到各个面的距离相等，是从六条棱出发的六个二面角的平分面的交点。且设四面体的表面积为S,内切球半径为R，则四面体的 体积为。 五、正三角形棱长为a时，周长为3a,面积为，高为，外接圆半径为，内切圆半径为。外接圆半径是内切圆半径的2倍。 正四面体棱长为a时，表面积为，高为，外接球半径为， 内切接球半径为。外接球半径是内切球半径的3倍。 六、任意三角形的三条中线交于一点，称为三角形的重心，重心到顶点的距离是它到对边中点距离的2倍。（重心定理）如图1所示：G为的重心。且 任意四面体的顶点与对面重心的连线交于一点，正是四面体的物理重心，且四面体的重心到顶点的距离是它到对面重心距离的3倍。（重心定理的推广） 如图2所示：E，F分别为的重心，AE与BF相交于点G，则G为四面体A-BCD的重心。 七、三角形中三个顶点的坐标分别为，

基于Matlab的相关频谱分析程序教程

Matlab 信号处理工具箱 谱估计专题 频谱分析 Spectral estimation （谱估计）的目标是基于一个有限的数据集合描述一个信号的功率（在频率上的）分布。功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信号的检测。 从数学上看，一个平稳随机过程n x 的power spectrum （功率谱）和correlation sequence （相关序列）通过discrete-time Fourier transform （离散时间傅立叶变换）构成联系。从normalized frequency （归一化角频率）角度看，有下式 ()()j m xx xx m S R m e ωω∞ -=-∞ = ∑ 注：()() 2 xx S X ωω=，其中()/2 /2 1 lim N j n n N n N X x e N ωω→∞=-=∑ πωπ-<≤。其matlab 近似为X=fft(x,N)/sqrt(N)，在下文中()L X f 就是指matlab fft 函数的计算结果了 使用关系2/s f f ωπ=可以写成物理频率f 的函数，其中s f 是采样频率 ()()2/s jfm f xx xx m S f R m e π∞ -=-∞ = ∑ 相关序列可以从功率谱用IDFT 变换求得： ()()()/2 2//2 2s s s f jfm f j m xx xx xx s f S e S f e R m d df f πωπ π ωωπ--= =? ? 序列n x 在整个Nyquist 间隔上的平均功率可以表示为 ()()() /2 /2 02s s f xx xx xx s f S S f R d df f π π ωωπ- -= =?? 上式中的

平面三角形单元有限元程序设计

平面三角形单元有限元程序设计

P 9 m 9 m 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=0.25，t=0.1m，忽略自重。试计算薄板的位移及应力分布。 要求： 1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3.提交程序编写过程的详细报告及计算机程序； 4.所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。

2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则 X Y P X Y P

边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎这样，由该方程求得的此位移值一定等于已知量。 （2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。 特点： （1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2）但这种方法不改变方程阶数，利于存贮。 （3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。 程序如下： 变量说明 NNODE 单元节点数 NPION 总结点数 NELEM 单元数 NVFIX 受约束边界点数 FIXED 约束信息数组 NFORCE 节点力数 FORCE 节点力数组

认识平面图形之三角形

认识平面图形之三角形 【知识框架】 1、图形分类（按不同标准给已知图形进行分类） 三角形的分类（认识直角三角形、锐角三角形、钝角三角形、等腰三角形、等边三角形） 2、三角形三角形内角和 3、三角形三边之间的关系 【知识要点】 图形分类三角形分类 1、把三角形按照不同的标准分类，并说明分类依据。 （1）按角分，分为：直角三角形、锐角三角形、钝角三角形，并了解其本质特征：三个角都是锐角的三角形是锐角三角形，有一个角是直角的三角形是直角三角形，有一个角是钝角的三角形是钝角三角形。 （2）按边分，分为：等腰三角形、等边三角形、任意三角形。有两条边相等的三角形是等腰三角形，三条边都相等的三角形是等边三角形。 2、通过分类，使学生弄清等腰三角形和等边三角形的关系：等边三角形是特殊 的等腰三角形。 三角形内角和 1、任意一个三角形内角和等于180度。 2、能应用三角形内角和的性质解决一些简单的问题。 三角形边的关系 1、三角形任意两边之和大于第三边。 2、根据上述知识点判断所给的已知长度的三条线段能否围成三角形。如果能围 成三角形，能围成一个什么样的三角形。 【公式概念】 1、围成三角形的条件：较短两条边长度的和一定大于第三条边。 2、从三角形的一个顶点到对边的垂直线段是三角形的高，这条对边是三角形的底。 3、三角形具有稳定性（也就是当一个三角形的三条边的长度确定后，这个三角形的形状和大小都不会改变），生活中很多物体利用了这样的特性。如：人字梁、斜拉桥、自行车车架。 4、三个角都是锐角的三角形是锐角三角形。（两个内角的和大于第三个内角。） 5、有一个角是直角的三角形是直角三角形。（两个内角的和等于第三个内角。两个锐角的和是90度。两条直角边互为底和高。） 6、有一个角是钝角的三角形是钝角三角形。（两个内角的和小于第三个内角。） 7、任意一个三角形至少有两个锐角，都有三条高，三角形的内角和都是180度。（锐角三角形的三条高都在三角形内；直角三角形有两条高落在两条直角边上；钝角三角形有两条高在三角形外）。 8、把一个三角形分成两个直角三角形就是画它的高。 9、两条边相等的三角形是等腰三角形，相等的两条边叫做腰，另外一条边叫做底，两条腰的夹角叫做顶角，底和腰的两个夹角叫做底角，它的两个底角也相等，是轴对称图形，有一条对称轴（跟底边高正好重合。）三条边都 第1页 相等的三角形是等边三角形，三条边都相等，三个角也都 相等（每个角都是60°，所有等边三角形的三个角都是60°。）

Matlab-PDE工具箱有限元法求解偏微分方程

在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。 偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。 偏微分方程基本形式 而以上的偏微分方程都能利用PDE工具箱求解。 PDE工具箱 PDE工具箱的使用步骤体现了有限元法求解问题的基本思路，包括如下基本步骤： 1) 建立几何模型 2) 定义边界条件 3) 定义PDE类型和PDE系数 4) 三角形网格划分

5) 有限元求解 6) 解的图形表达 以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上，如下 具体实现如下。 打开工具箱 输入pdetool可以打开偏微分方程求解工具箱，如下 首先需要选择应用模式，工具箱根据实际问题的不同提供了很多应用模式，用户可以基于适

当的模式进行建模和分析。 在Options菜单的Application菜单项下可以做选择，如下 或者直接在工具栏上选择，如下 列表框中各应用模式的意义为： ① Generic Scalar：一般标量模式（为默认选项）。 ② Generic System：一般系统模式。 ③ Structural Mech.，Plane Stress：结构力学平面应力。 ④ Structural Mech.，Plane Strain：结构力学平面应变。

⑤ Electrostatics：静电学。 ⑥ Magnetostatics：电磁学。 ⑦ Ac Power Electromagnetics：交流电电磁学。 ⑧ Conductive Media DC：直流导电介质。 ⑨ Heat Tranfer：热传导。 ⑩ Diffusion：扩散。 可以根据自己的具体问题做相应的选择，这里要求解偏微分方程，故使用默认值。此外，对于其他具体的工程应用模式，此工具箱已经发展到了Comsol Multiphysics软件，它提供了更强大的建模、求解功能。 另外，可以在菜单Options下做一些全局的设置，如下 l Grid：显示网格 l Grid Spacing…：控制网格的显示位置 l Snap：建模时捕捉网格节点，建模时可以打开 l Axes Limits…：设置坐标系围 l Axes Equal：同Matlab的命令axes equal命令 建立几何模型 使用菜单Draw的命令或使用工具箱命令可以实现简单几何模型的建立，如下 各项代表的意义分别为

平面图形的认识---三角形的认识综合提优(压轴题)

平面图形的认识---三角形的认识综合提优 如图，在△ABC中，AD平分∠BAC，P为线段AD上的一个动点，PE⊥AD交直线BC于点 E. (1)若∠B=35°,∠ACB=85°，求∠E的度数； (2)当P点在线段AD上运动时，猜想∠E与∠B、∠ACB的数量关系，写出结论无需证明。 已知如图，射线CB∥OA，∠C=∠OAB=100°，E、F在CB上，且满足∠FOB=∠AOB，OE 平分∠COF。 （1）求∠EOB的度数； （2）若平行移动AB，那么∠OBC∶∠OFC的值是否随之变化？若变化，找出变化规律；若不变，求出这个比值； （3）在平行移动AB的过程中，是否存在某种情况，使∠OEC=∠OBA？若存在，求出其度数；若不存在，说明理由。 如图，AD⊥BD，AE平分∠BAC，∠B＝30°，∠ACE＝110°．求∠AED的度数．

现有两块大小相同的直角三角板△ABC、△DEF，∠ACB＝∠DFE＝90°，∠A＝∠D＝30°． (1)将这两块三角板摆成如图①的形式，使B、F、E、A在同一条直线上，点C在边DF上，DE与AC相交于点G，试求∠AGD的度数； (2)将图①中的△ABC固定，把△DEF绕着点F逆时针旋转成如图②的形式，当旋转 的角度等于多少度时，DF∥AC?并说明理由． 如图，△ABC中，∠ABC=50°，∠ACB=70°，D为边BC上一点（D与B、C不重合），连接AD，∠ADB的平分线所在直线分别交直线AB、AC于点E、F． （1）求证：2∠AED-∠CAD=170°； （2）若∠ABC=∠ACB=n°，且D为射线CB上一点，（1）中其他条件不变，请直接写出∠AED与∠CAD的数量关系．（用含n的代数式表示）