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

. 一、题目

如图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

适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。

做法：

（1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同

时将载荷列阵｛R｝中相应元素用已知位移置换。

◎这样，由该方程求得的此位移值一定等于已知量。

（2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。

◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为

保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。

特点：

（1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。

（2）但这种方法不改变方程阶数，利于存贮。

（3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。

程序如下：

变量说明

NNODE 单元节点数

NPION 总结点数

NELEM 单元数

NVFIX 受约束边界点数

FIXED 约束信息数组

NFORCE 节点力数

FORCE 节点力数组

COORD 结构节点坐标数组

LNODS 单元定义数组

YOUNG 弹性模量

POISS 泊松比

THICK 厚度

B 单元应变矩阵(3*6)

D 单元弹性矩阵(3*3)

S 单元应力矩阵(3*6)

A 单元面积

ESTIF 单元刚度矩阵

ASTIF 总体刚度矩阵

ASLOD 总体荷载向量

ASDISP 节点位移向量

ELEDISP 单元节点位移向量

STRESS 单元应力

%********************************************************** %初始化

clear

format short e %设定输出类型

clear %清除内存变量

NELEM=36 %单元个数（单元编码总数）

NPION=28 %结点个数（结点编码总数）

NVFIX=2 %受约束边界点数

NFORCE=1 %结点荷载个数

YOUNG=2e11 %弹性模量

POISS=0.25 %泊松比

THICK=0.1 %厚度

LNODS=[1 2 3;2 4 5;2 5 3;3 5 6;

4 7 8;4 8 5;

5 8 9;5 9 6;

6 9 10;

7 11 12;7 12 8;

8 12 13;

8 13 9;9 13 14;9 14 10;10 14 15;

11 16 17;11 17 12; 12 17 18; 12 18 13;

13 18 19; 13 19 14;14 19 20;14 20 15;

15 20 21;16 22 23;16 23 17;17 23 24;

17 24 18;18 24 25;18 25 19;25 19 26;

19 26 20;20 26 27;20 27 21;21 27 28] %单元定义数组（单元结点号）

%相应为单元结点号（编码）、按逆时针顺序输入

COORD=[0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3;

1.5 3;-

2.25 4.5;-0.75 4.5;0.75 4.5;

2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6;

-3.75 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5;

2.25 7.5;

3.75 7.5;-

4.5 9;-3 9;

-1.5 9;0 9;1.5 9;3 9;4.5 9] %结点坐标数组

%坐标：x，y 坐标（共NPOIN 组）

FORCE=[1 0 -15] %结点力数组（受力结点编号, x 方向,y 方向）FIXED=[22 1 1;28 1 1] %约束信息（约束点，x 约束，y 约束）

%有约束为1，无约束为0

%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0

%********************************************************** for i=1:NELEM

%生成弹性矩阵 D

D= [1 POISS 0;

POISS 1 0;

0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)

%********************************************************** %计算当前单元的面积

A=-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

%********************************************************** %生成应变矩阵 B

for j=0:2

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

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

end

B=[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*A);

B1( :,:,i)=B;

%********************************************************** %求应力矩阵S=D*B

S=D*B;

ESTIF=B'*S*THICK*A; %求解单元刚度矩阵

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

for j=1:3

for k=1:3

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF((a(j)*2-1) :a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);

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

%度矩阵中

end

end

end

%********************************************************** %将约束信息加入总体刚度矩阵（对角元素改一法）

for i=1:NVFIX

if FIXED(i,2)==1

ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零

ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零

ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; %对角元素为 1

end

%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵

%**********************************************************

if FIXED(i,3)==1

ASTIF( :,FIXED(i,1)*2)=0; %一列为零

ASTIF(FIXED(i,1)*2,:)=0; %一行为零

ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为 1

end

end

%********************************************************** %生成荷载向量

ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCE

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

end

%********************************************************** %求解内力

ASDISP=ASTIF\ASLOD' %计算节点位移向量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

i

STRESS=D*B1(:, :, i)*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

有限元2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)综述

第2章 弹性力学平面问题有限单元法 2.1 三角形单元(triangular Element) 三角形单元是有限元分析中的常见单元形式之一，它的优点是: ①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。 一、结点位移和结点力列阵 设右图为从某一结构中取出的一典型三角形单元。 在平面应力问题中，单元的每个结点上有沿x 、y 两个方向的力和位移，单元的结点位移列阵规定为： 相应结点力列阵为: (式2-1-1) 二、单元位移函数和形状函数 前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构 造）一组在单元内有定义的位移函数作为近似计算的基础。即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。 构造位移函数的方法是：以结点(i,j,m)为定点。以位移(u i ,v i ,…u m v m )为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。 在平面应力问题中，有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成: (,)123 u u x y x y ααα==++ 546(,)v v x y x y ααα==++ (2-1-2)a 式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标) {}??? ?? ?????=????? ???? ?????????????=m j i m e d d d d m j j i v u v u v u i {} i i j j m X Y X (2-1-1)Y X Y i e j m m F F F F ?? ?? ???? ???? ??==??????????????????

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

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

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

第三章平面问题的有限元法作业及答案

第三章 平面问题的有限元法作业 1. 图示一个等腰三角形单元及其节点编码情况，设μ=0，单元厚度为t 。求 1）形函数矩阵[]N ；2）应变矩阵[]B ；3）应力矩阵[]S 。 4 第1题图 第2题图 2. 如题图所示，结构为边长等于a 的正方形，已知其节点位移分别为：11(,)u v 、 22(,)u v 、33(,)u v 、44(,)u v 。试求A 、B 、C 三点的位移。其中A 为正方形形心，B 为三角形形心。 3.直角边边长为l 的三角形单元，如题图所示。试计算单元等效节点载荷列阵（单元厚度为t ，不计自重）。 第3题图 第4题图 4. 如题图所示，各单元均为直角边边长等于l 的直角三角形。试计算（1）单元等效节点载荷列阵；（2）整体等效节点载荷列阵。已知单元厚度为t ，不计自重。

5．下列3个有限元模型网格，哪种节点编号更合理？为什么？ 9 34 6 7912 11 34 6 12142 (a) (b) (c) 第5题图 6.将图示结构画出有限元模型；标出单元号和节点号；给出位移边界条件；并计算半带宽（结构厚度为t ）。 2a (a) (b) 无限长圆筒 (c) 第6题图 7. 结构如图所示，已知结构材料常数E 和 ，单元厚度为t 。利用结构的对称性，采用一个单元，分别计算节点位移和单元应力。 第7题图

答案： 1. 1）形函数 i x N a = ， j y N a = ， 1m x y N a a =-- 2）应变矩阵 []1000101 000101011011B a -????=-??--???? 3）应力矩阵 []100010100 01 0111 110022 2 2S a ? ???-? ?=-????- -? ?? ? 2. A 点的位移为 ()2312A u u u = + ， ()231 2A v v v =+ B 点的位移为 ()24313B u u u u = ++ ， ()2431 3B v v v v =++ C 点的位移为 ()1223C a u u u = + ， ()C 1223 a v v v =+ 3. 单元等效节点载荷列阵为 {}11 11 00003 663 T e i j i j R q q q q ?? =++?? ?? 4. （2）整体等效节点载荷向量为 {}111100006 322T R qlt P qlt P P qlt qlt ?? =-???? 7. （1） 减缩后的整体刚度方程 22 12 2 1222 22221110222021102(1)2 2102x x b b ab R b ab b P v Et ab a b ab ab R v b a μμμ μμμμμμ---??- - ??????????--?????? -??? ?=????---+ +? ???? ?????????-????+?? ? ? 节点位移

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

平面三角形单元有限元程序设计 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 节点编号 单元编号

平面问题的有限元法-Read

3 弹性力学平面问题的有限元法 本章包括以下的内容： 3.1弹性力学平面问题的基本方程 3.2单元位移函数 3.3单元载荷移置 3.4单元刚度矩阵 3.5单元刚度矩阵的性质与物理意义 3.6整体分析 3.7约束条件的处理 3.8整体刚度矩阵的特点与存储方法 3.9方程组解法 3.1弹性力学平面问题的基本方程 弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程，把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。 弹性力学的基本假定如下： 1）完全弹性，2）连续，3）均匀，4）各向同性，5）小变形。 3.1.1基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变，各自的定义如下。 体力 体力是分布在物体体积内的力，例如重力和惯性力。 面力 面力是分布在物体表面上的力，例如接触压力、流体压力。 应力 物体受到约束和外力作用，其内部将产生内力。物体内某一点的内力就是应力。 图3.1

如图3.1假想用通过物体内任意一点p 的一个截面mn 将物理分为Ⅰ、Ⅱ两部分。将部分Ⅱ撇开，根据力的平衡原则，部分Ⅱ将在截面mn 上作用一定的内力。在mn 截面上取包含p 点的微小面积A ?，作用于A ?面积上的内力为Q ?。 令A ?无限减小而趋于p 点时，Q ?的极限S 就是物体在p 点的应力。 S A Q A =??→?0lim 应力S 在其作用截面上的法向分量称为正应力，用σ表示；在作用截面上的切向分量称为剪应力，用τ表示。 显然，点p 在不同截面上的应力是不同的。为分析点p 的应力状态，即通过p 点的各个截面上的应力的大小和方向，在p 点取出的一个平行六面体，六面体的各楞边平行于坐标轴。 图3.2 将每个上的应力分解为一个正应力和两个剪应力，分别与三个坐标轴平行。用六面体表面的应力分量来表示p 点的应力状态。应力分量的下标约定如下： 第一个下标表示应力的作用面，第二个下标表示应力的作用方向。 xy τ，第一个下标x 表示剪应力作用在垂直于X 轴的面上，第二个下标y 表示剪应力指 向Y 轴方向。 正应力由于作用表面与作用方向垂直，用一个下标。x σ表示正应力作用于垂直于X 轴的面上，指向X 轴方向。 应力分量的方向定义如下： 如果某截面上的外法线是沿坐标轴的正方向，这个截面上的应力分量以沿坐标轴正方向为正； 如果某截面上的外法线是沿坐标轴的负方向，这个截面上的应力分量以沿坐标轴负方向为正。 剪应力互等：xz zx zy yz yx xy ττττττ===,, 物体内任意一点的应力状态可以用六个独立的应力分量x σ、y σ、z σ、xy τ、yz τ、zx τ

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

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

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 节点力数组

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

三角形单元有限元程序设计 %******************************************************************* % NNODE 单元节点数NPION 总结点数 % NELEM 单元数NVFIX 受约束边界点数 % FIXED 约束信息数组NFORCE 节点力数 % FORCE 节点力数组COORD 结构节点坐标数组 % LNODS 单元定义数组 % YOUNG 弹性模量POISS 泊松比 % THICK 厚度 B 单元应变矩阵(3*6) % D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6) % A 单元面积ESTIF 单元刚度矩阵 % ASTIF 总体刚度矩阵ASLOD 总体荷载向量 % ASDISP 节点位移向量ELEDISP 单元节点位移向量 % STRESS 单元应力FP1 数据文件指针 %******************************************************************* format short e %设定输出类型 clear %清除内存变量 FP1=fopen('1-1.txt','rt'); %打开输入数据文件存放初始数据 %读入控制数据 NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数） NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数） NVFIX=fscanf(FP1,'%d',1) %受约束边界点数 NFORCE=fscanf(FP1,'%d',1), %结点荷载个数 YOUNG=fscanf(FP1,'%e',1), %弹性模量 POISS=fscanf(FP1,'%f',1), %泊松比 THICK=fscanf(FP1,'%f',1) %厚度 LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号） %相应为单元结点号（编码）、按逆时针顺序输入 COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组 %坐标：x，y坐标（共NPOIN组） FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组 %(n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向 FIXED=fscanf(FP1,'%d',[3,NVFIX])' %约束信息数组 % (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0 %******************************************************************* %生成单元刚度矩阵并组成总体刚度矩阵 ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 %*******************************************************************

弹性力学平面问题的有限元法

Mmm 3 弹性力学平面问题的有限元法 本章包括以下的内容： 3.1弹性力学平面问题的基本方程 3.2单元位移函数 3.3单元载荷移置 3.4单元刚度矩阵 3.5单元刚度矩阵的性质与物理意义 3.6整体分析 3.7约束条件的处理 3.8整体刚度矩阵的特点与存储方法 3.9方程组解法 3.1弹性力学平面问题的基本方程 弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程，把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。 弹性力学的基本假定如下： 1）完全弹性，2）连续，3）均匀，4）各向同性，5）小变形。 3.1.1基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变，各自的定义如下。 体力 体力是分布在物体体积内的力，例如重力和惯性力。 面力 面力是分布在物体表面上的力，例如接触压力、流体压力。 应力 物体受到约束和外力作用，其内部将产生内力。物体内某一点的内力就是应力。

图3.1 如图3.1假想用通过物体内任意一点p 的一个截面mn 将物理分为Ⅰ、Ⅱ两部分。将部分Ⅱ撇开，根据力的平衡原则，部分Ⅱ将在截面mn 上作用一定的内力。在mn 截面上取包含p 点的微小面积A ?，作用于A ?面积上的内力为Q ?。 令A ?无限减小而趋于p 点时，Q ?的极限S 就是物体在p 点的应力。 S A Q A =??→?0lim 应力S 在其作用截面上的法向分量称为正应力，用σ表示；在作用截面上的切向分量称为剪应力，用τ表示。 显然，点p 在不同截面上的应力是不同的。为分析点p 的应力状态，即通过p 点的各个截面上的应力的大小和方向，在p 点取出的一个平行六面体，六面体的各楞边平行于坐标轴。 图3.2 将每个上的应力分解为一个正应力和两个剪应力，分别与三个坐标轴平行。用六面体表面的应力分量来表示p 点的应力状态。应力分量的下标约定如下： 第一个下标表示应力的作用面，第二个下标表示应力的作用方向。 xy τ，第一个下标x 表示剪应力作用在垂直于X 轴的面上，第二个下标y 表示剪应力指 向Y 轴方向。

有限元作业：三角形单元求解

《有限元作业》 年级2015级 学院机电工程学院 专业名称 班级学号 学生姓名 2016年05月

如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个） 图(一) 图（二）三角形单元求解 图（三）四边形单元求解

（1）如图划分三角形单元，工分成四个分别为???④ （2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N 调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1 = +06 * 0 0 0 0 0 0 0 0 k2 = +06 * 0 0 0 0 0 0 0 0 k3 = +06 * 0 0 0 0 0 0 0 0 k4 = +06 *

0 0 0 0 0 0 0 0 调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵 求出的节点位移 U = 调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy S1 = +03 * S2 = +03 *

S3 = +03 * S4 = +03 * 二、 （1）如图划分四边形单元，工分成四个分别为?? （2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N 调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵 调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵

有限元三角形等参单元

北方工业大学 高等有限元课程总结 姓名:韩双鹏学号: 2011303310105 专业班级:结构研-11 系(部、院):建筑工程学院 2012 年5 月25 日

高等有限元学习总结——六节点三角形等参数单元 1 概述 从弹性力学基本方程到有限元原理再到最新进展，经过本课程的学习，比较系统的掌握了有限元相关内容，更学习到了一种方法、一些生活中的哲理。首先从大方向掌握所学内容，避免迷失在局部造成一叶遮目不见泰山之悲剧，比如弹性力学原理从大方向说就是三类方程，以及其在各类问题中的应用；其次了解了科研的相关过程及创新之处，从已知的东西到无知的领域，正如老师所说，能成功地把某一领域的东西搬到相关领域，这就是一大创造，比如有限元中将梁弯曲的理论研究厚板弯曲问题，由有限元标准单元到等参元的研究等；再有，我们生活中的常识、学习中的某些东西值得我们细细品味，也许这就是平时所说的小事反应大道理，老师的理论：“很多想法都是错误的”“很好想到的方法也许很难走通”“有缺陷的东西才更体现出美”“平衡的理论，吃点亏也许是福”等等，受益匪浅。不再一一赘述，本文将取其中的一个知识点，总结六节点三角形等参单元的相关内容。 我们知道，无论三节点或者六节点三角形单元还是四节点或者八节点矩形单元，它们形状简单、规则但计算精度低，且对于复杂边界的适应性差，难以很好的拟合曲边边界，解决这一问题的通用方法是细分边界，以直代曲，利用更多的简单单元去拟合边界复杂的区域。但这样处理仍存在折线代替曲线所带来的误差，且这种误差不能通过提高单元位移函数的精度来补偿。那么能否构造出单元形状任意、边界适应性好、计算精度高的曲边单元，以便在给定的精度下用较少数目的单元去解决实际问题？这就是有限元中一类重要的单元——等参数单元。本文将总结等参数单元的基本概念，并以六节点三角形单元为例讲述等参元实现过程中的三种变换，以及该等参元的收敛性等问题。 2 等参数单元及实现过程 2.1 等参数单元概念 由于实际问题的复杂性，通常需要使用一些形状不规整和形状复杂的单元来离散边界形状复杂的原问题。如下图所示（a）中为常见的几何形状不规整的实际单元，称为实际单元，也称为参数单元。（b）中为对应的形状规整的单元，称为标准单元。对于形状复杂的实际单元的单元分析，若仍采用前面介绍的方法进行，则在单元位移函数的建立和单元刚度矩阵计算方面会遇到许多困难。由此可考虑利用前面介绍过的形状规整的标准单元的单元分析来研究实际单元，几何形状的不同可认为是坐标变换的结果。

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

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

网格划分 如图，将薄板如图划分为6行，并建立坐标系，则 X Y P X Y P 节点编号 单元编号

刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第奇/偶数个计算 单元刚度的集成： [][] [][] [][] ' ' ' ' ' ' 2 1 56 56 6 6 56 56 2 6 6 2 56 56 1 6 6 1 e Z e e e Z e Z e e e e k k k K k k k k k k + ? + + = ? = ? = = ? = = ? = ? ? ? ? ? ? 边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎这样，由该方程求得的此位移值一定等于已知量。 （2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。 特点： （1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2）但这种方法不改变方程阶数，利于存贮。

有限元平面问题

[] ()1 ,2 1 1112 1 =∴---++== i i i j k i j k i i k k j j i k k j j i i y x N y x y x y x y x y x y x y x y x y x A 即：()()()1,,,===k k k j j j i i i y x N y x N y x N （由i,j,k 轮换性知） 同理可证：()()0,,==k k i j j i y x N y x N （作业：证明：()()k j i j i y x N j j i ,,0,=≠=） 因此 ()()()()()()()()()()?? ??? ??===???≠===?======1 ,,0,,0,,1,00,,1,,0,0 ,,0,,1,k k k j j k i i k j i i k k j j j j i i j k k i j j i i i i y x N y x N y x N j i j i j N y x N y x N y x N y x N y x N y x N δ （2-12） 即形函数在自己节点上为1，在其余节点上为0。 2. 在单元上任意一点，三个形函数之和为1，即 ()()()1,,,=++y x N y x N y x N k j i 。 证明： ()()()()()()()[] y c c c x b b b a a a A y c x b a y c x b a y c x b a A y x N y x N y x N k j i k j i k j i k k k j j j i i i k j i ++++++++= ++++++++=++21 21 ,,,

有限元作业：三角形单元求解

有限元作业：三角形单元求 解 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

《有限元作业》 年级2015级 学院机电工程学院 专业名称 班级学号 学生姓名 2016年05月

如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比 V=1/6,厚度为I=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个） 图(一) 图（二）三角形单元求解

图（三）四边形单元求解（1）如图划分三角形单元，工分成四个分别为 ④（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N 调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1 = 1.0e+06 * 7.2857 -3.0000 -2.1429 0.8571 -5.1429 2.1429 -3.0000 7.2857 2.1429 -5.1429 0.8571 -2.1429 -2.1429 2.1429 2.1429 0 0 -2.1429 0.8571 -5.1429 0 5.1429 -0.8571 0 -5.1429 0.8571 0 -0.8571 5.1429 0 2.1429 -2.1429 -2.1429 0 0 2.1429 k2 = 1.0e+06 * 5.1429 0 -5.1429 0.8571 0 -0.8571 0 2.1429 2.1429 -2.1429 -2.1429 0 -5.1429 2.1429 7.2857 -3.0000 -2.1429 0.8571 0.8571 -2.1429 -3.0000 7.2857 2.1429 -5.1429 0 -2.1429 -2.1429 2.1429 2.1429 0

平面问题的有限元方法-8.6作业-集中力-平面问题

一个受到集中力P 作用的结构，泊松比ν=61，m N P y 16=，t=1cm ，试按平面应力问题计算，采用三角形单元，求出节点位移。 解： 如图所示，划分三角形单元为四部分，并进行单元坐标编号，编程进行求解

单元①的刚度矩阵为： ???? ??????=333231232221 131211 1K K K K K K K K K K ()3,2,1===m j i 其中子矩阵表达式为： ???? ??????-+-+-+-+?-=s r s r s r s r s r s r s r s r rs b b v c c c b b c b c c b c c v b b Et K 21212121)1(42ννννν()m j i s r ,,,= E E E E E Et 2,,22210944.110944.121)611(401 .0)1(4--==?≈???? ??-??=? -ν 调用 Triangle2D3Node_Stiffness 函数，求出单元刚度矩阵： )3,2,1(0.2143 0 0 0.2143- 0.2143- 0.2143 0 0.5143 0.0857- 0 0.0857 0.5143- 0 0.0857- 0.5143 0 0.5143- 0.0857 0.2143- 0 0 0.2143 0.2143 0.2143- 0.2143- 0.0857 0.5143- 0.2143 0.7286 0.3000- 0.2143 0.5143- 0.0857 0.2143- 0.3000- 0.7286 '1===????????? ???????????=m j i E K

有限元平面问题MATLAB程序

计算力学(有限元方法部分) 程序设计 程序说明书 程序1：平面问题的有限元分析 文件名：h01.m 算例文本：h01.txt 输出文本：result1.txt 使用前请先将h01.txt放入默认的文本读取路径（我的要求与m文件在同一文件夹内）！ 文本输入顺序： 材料信息（编号、弹性模量、泊松比） 注意：材料编号须按1、2、3、4……的顺序排列 节点信息（编号、x坐标、y坐标） 注意：节点编号须按1、2、3、4……的顺序排列 约束信息（约束节点号、x方向有无约束、y方向有无约束、x方向位移、y 方向位移）有约束处填一正数，无约束处填0。无约束处请勿输入位移。 单元信息（厚度、材料编号、节点编号，若为3节点单元，则第四个编号填0） 注意：单元编号须按1、2、3、4……的顺序排列 集中力（作用节点号、x方向力、y方向力） 线荷载（作用边上的两个节点号、x方向分布力、y方向分布力） 面荷载（作用单元号、x方向分布力、y方向分布力） 程序可调部分： 4-6行中可以指定输出哪些图像（按顺序依次为节点、单元图像，x、y方向位移图像，xx、yy、xy方向应力图像），第7行中可以指定输入的.txt文本名称。 文本输出内容： 结点位移信息（节点号、x方向位移、y方向位移） 单元形心处的应变信息（单元号、x方向正应变、y方向正应变、xy方向工程切应变）

单元形心处的应力信息（单元号、x方向正应力、y方向正应力、xy方向切应力） 本程序附有三角形单元自动加密前处理部分h01auto.m，其算例文本： h01coarse.txt，输出文本：h01new.txt。它可以适用于本题的要求，在已给定本题所需全部信息的情况下将已有的单元加密为三角形单元。其输出文本可直接作为上面程序的输入文本。 h01.m h01.txt h01auto.m h01coarse.txt 欢迎交流与提问！附上邮箱：x67891@https://www.sodocs.net/doc/023539023.html,。祝力学学习顺利！

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

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

刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算

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

、题目 如图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行，并建立坐标系，则 /2\/4\ 10xl0 八18八、20 人？2/\24 「八刘人却八引八読人悶 /、

节点编号 单兀编号 刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1 ）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第奇/偶数个计算 k1e66 k1 ' 56 56 k2e66 k2 ' 56 56 单元刚度的集成：kZ66 k Z '56 56 K k e' k e'12 k Z e' 边界约束的处理：划0 置1 法 适用：这种方法适用于边界节点位移分量为已知(含为0) 的各种约束。 做法：

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

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

刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算

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

一、题目 如图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