常应变三角形单元

合集下载

三角形常应变单元程序的编制与使用11页

三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。

有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。

对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。

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 THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1●说明：NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y 方向；FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)，DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。

有限元基础知识归纳

有限元知识点归纳1.、有限元解的特点、原因？答：有限元解一般偏小，即位移解下限性原因：单元原是连续体的一局部，具有无限多个自由度。

在假定了单元的位移函数后，自由度限制为只有以节点位移表示的有限自由度，即位移函数对单元的变形进行了约束和限制，使单元的刚度较实际连续体加强了，因此，连续体的整体刚度随之增加，离散后的刚度较实际的刚度K为大，因此求得的位移近似解总体上将小于精确解。

2、形函数收敛准那么〔写出某种单元的形函数，并讨论收敛性〕P49(1)在节点i处N i=1,其它节点N i=0;(2)在单元之间，必须使由其定义的未知量连续；(3)应包含完全一次多项式；(4)应满足∑Ni=1以上条件是使单元满足收敛条件所必须得。

可以推证，由满足以上条件的形函数所建单元是完备协调的单元，所以一定是收敛的。

4、等参元的概念、特点、用时注意什么？〔王勖成P131〕答：等参元—为了将局部坐标中几何形状规那么的单元转换成总体〔笛卡尔〕坐标中的几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立一个坐标变换。

即：为建立上述的变换，最方便的方法是将上式表示成插值函数的形式，即：其中m是用以进行坐标变换的单元节点数，xi,yi,zi是这些结点在总体〔笛卡尔〕坐标内的坐标值，Ni’称为形状函数，实际上它也是局部坐标表示的插值函数。

称前者为母单元，后者为子单元。

还可以看到坐标变换关系式和函数插值表示式：在形式上是相同的。

如果坐标变换和函数插值采用相同的结点，并且采用相同的插值函数，即m=n，Ni’=Ni,那么称这种变换为等参变换。

5、单元离散？P42答：离散化既是将连续体用假想的线或面分割成有限个局部，各局部之间用有限个点相连。

每个局部称为一个单元，连接点称为结点。

对于平面问题，最简单、最常用的离散方式是将其分解成有限个三角形单元，单元之间在三角形顶点上相连。

这种单元称为常应变三角形单元。

常用的单元离散有三节点三角形单元、六节点三角形单元、四节点四边形单元、八节点四边形单元以及等参元。

计算力学(有限单元法)第三章重点整理

第三章一、三角形单元（常应变单元）1）三角形单元位移函数：123456u a a x a yv a a x a y =++⎧⎨=++⎩2）位移函数用形函数来表示：i i j j m mi i j j m mu N u N u N u v N v N v N v =++⎧⎨=++⎩其中1()(,,)2i i i i N a bx c y i j m A =++，,(,,)i j m m ji j m ij m a x y x y b y y i j m c x x ⎧=-⎪=-⎨⎪=-+⎩，11121i i j j mmx y A x y x y = 形函数用单元节点位移分量来描述位移函数的插值函数，反映了单元的位移形态，数学是反映了节点位移对单元内任一点位移的插值。

矩阵形式：0000i i ijm j ijm jm m u v N N N u u N NN v v u v ⎧⎫⎪⎪⎪⎪⎪⎪⎡⎤⎧⎫⎪⎪=⎨⎬⎨⎬⎢⎥⎩⎭⎣⎦⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭或{}[]{}[][]{}i j m f N N N N δδ⎡⎤⎡⎤==⎣⎦⎣⎦4）单元应变：{}[]{}B εδ=（其中[]B 为常量）由x y xy u x v y u v y x εεγ⎧⎫∂⎪⎪∂⎪⎪⎧⎫⎪⎪∂⎪⎪=⎨⎬⎨⎬∂⎪⎪⎪⎪⎩⎭⎪⎪∂∂+⎪⎪∂∂⎩⎭得到[]001002ii i i i i i i i Nx b N B c y Ac b N N yx ⎡⎤∂⎢⎥∂⎢⎥⎡⎤⎢⎥∂⎢⎥==⎢⎥⎢⎥∂⎢⎥⎢⎥⎣⎦⎢⎥∂∂⎢⎥∂∂⎣⎦应变和节点位移关系式：00010002i i x i j m j y i j m j xy iijjmm m m u v b b b u c c c v A c b c b c b u v εεγ⎧⎫⎪⎪⎪⎪⎧⎫⎡⎤⎪⎪⎪⎪⎢⎥⎪⎪=⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪⎢⎥⎩⎭⎣⎦⎪⎪⎪⎪⎪⎪⎩⎭5）单元应力：{}[][]{}{}[]D B S σδδ==其中36[][[][][]]i j k S S S S ⨯=平面应力问题2[],(,,)2(1)1122i i i i ii i b c ES b c i j m Ac b μμμμμ⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥--⎢⎥⎣⎦平面应变问题将上式中的21E E μ=-6）单元平衡方程：{}{}[]d k F δ=,{}{}{}{}d V S c F F F p =++7）单元刚度矩阵：[][][][]TVk B D B dv=⎰（表示单元力和单元位移关系间的系数，代表单元的刚度特性）性质：（1）三角形单元刚度矩阵与坐标系无关，即单元刚度矩阵[]k 不随单元或坐标轴的平行移动或作n π角度的转动而改变（平面问题的单元刚度矩阵可以认为是结构坐标系中的单元刚度矩阵，没有坐标变换问题）（2）单元刚度矩阵中每个元素ij k 的物理意义表示单元第j 个自由度产生单位位移，其它自由度固定时，第i 个自由度产生的节点力。

4平面问题有限元分析

平面问题有限元分析
引言常应变三角形单元矩形双线性单元三角形类单元形函数矩形类单元形函数平面等参数单元 Wilson 非协调元及程序
引
言
杆系问题以结点作为分割单元的“结点”是很自然的，杆系问题以结点作为分割单元的“结点”是很自然的，但对于平面问题，待分析物体是连续的，但对于平面问题，待分析物体是连续的，并不存在实际结点。要将物体“ 成单元，结点。要将物体“拆”成单元，必须用一些假想的线或将物体进行分割时，将物体进行分割时，必须保证相面作人为地分割。面作人为地分割。邻单元具有公共边界。假定相邻单元仅在一些点（顶点邻单元具有公共边界。假定相邻单元仅在一些点（或顶点加边中点）相连接。这些点即为“结点” 或顶点加边中点）相连接。这些点即为“结点”。实际计算时，可将连续体分成多种形状单元，为讨论简单，计算时，可将连续体分成多种形状单元，为讨论简单，现暂时规定只用一种单元来分割。现暂时规定只用一种单元来分割。以位移为未知量的有限元法，以位移为未知量的有限元法，最关键的工作是建立单元位移场，因此本章主要介绍各种单元位移场的建立。元位移场，因此本章主要介绍各种单元位移场的建立。平面问题有限元法可用的单元很多，平面问题有限元法可用的单元很多，先介绍两种最简单的单元：三角形和矩形。然后再介绍其它的单元。单的单元：三角形和矩形。然后再介绍其它的单元。
常应变三角形单元
由于面积坐标有形函数性质， 3 位移模式由于面积坐标有形函数性质，因 3 此根据试凑法可得形函数= 形函数 Ni=Li = 面积坐标 y P 位移为u 如果结点 i 位移为 i、vi，则 2 单元位移模式（位移场）单元位移模式（位移场）为 1 x u=Σ Niui ； v=Σ Nivi Σ Σ 1) 面积坐标和直角坐标关系

结构力学练习题及答案

一．是非题（将判断结果填入括弧：以O 表示正确，X 表示错误）（本大题分4小题，共11分）1 . (本小题 3分)图示结构中DE 杆的轴力F NDE =F P /3。

（）.2 . (本小题 4分)用力法解超静定结构时，只能采用多余约束力作为基本未知量。

（）3 . (本小题 2分)力矩分配中的传递系数等于传递弯矩与分配弯矩之比，它与外因无关。

（）4 . (本小题 2分)用位移法解超静定结构时，基本结构超静定次数一定比原结构高。

（）二．选择题（将选中答案的字母填入括弧内）（本大题分5小题，共21分） 1 （本小题6分）图示结构EI=常数，截面A 右侧的弯矩为：（）A ．2/M ；B ．M ；C ．0； D. )2/(EI M 。

2. （本小题4分）图示桁架下弦承载，下面画出的杆件内力影响线，此杆件是：（）Ａ.ch; B.ci; C.dj; D.cj.F p /2M2a2a a aa aA F p /2F p /2 F p /2F p F pa a aa F PED3. (本小题 4分)图a 结构的最后弯矩图为：A. 图b;B. 图c;C. 图d;D.都不对。

（）( a) (b) (c) (d)4. (本小题 4分)用图乘法求位移的必要条件之一是： A.单位荷载下的弯矩图为一直线； B.结构可分为等截面直杆段； C.所有杆件EI 为常数且相同； D.结构必须是静定的。

( ) 5. （本小题3分）图示梁A 点的竖向位移为（向下为正）：( ) Ａ．F P l 3/(24EI); B. F P l 3/(!6EI); C. 5F P l 3/(96EI); D. 5F P l 3/(48EI).三（本大题 5分）对图示体系进行几何组成分析。

A l /2l /2EI 2EIF Pa d c eb fgh iklF P =11j llM /4 3M /4M /43M /43M /4M /4M /8 M /2EIEIM四（本大题 9分）图示结构B 支座下沉4 mm ，各杆EI=2.0×105 kN ·m 2，用力法计算并作M 图。

9第2章弹性力学平面问题及空间问题有限元

v u v 2 ， y 6 ， xy 3 5 都是常量，即线性位移模式反映 x y y x
假定的位移函数是多项式，它是连续函数，可以肯定，在单元内部位移函数是单值连续的。由于单元的位移函数 u 、 v 都是坐标 x 、 y 的线性函数，在单元边界上位移也是线性变化的，两个相邻单元在公共节点上具有相同的节点位移，因而相邻单元在公共边界上位移连续，即协调条件得到满足。由上面分析可以看出，三角形常应变单元的位移模式可以保证计算结果的收敛。
px
py
px
py ]
T
（2-1-7b）
（2 ）若在 jm 边上受线性分布的水平方向的面力，它在 j 点的集度为 q ，在 m 点的集度为零（如图 2-5）。可预计由该面力求得的等效节点载荷只有 R xj 、
R xm ，其余节点载荷分量必为零。
将 jm 边上的分布面力写成 s 的函数，为
s { p} [ (1 ) q 0]T l 在 jm 边上的形函数也需用变量 s 表示，根据形函数的含义，
Ve
[k ii ] [k ij ] [ k im ] [k ji ] [k ij ] [k jm ] [k mi ] [ k mj ] [k mm ]
式中， t 为单元的厚度，当单元划分得足够小时，可以认为每个单元的厚度 t 为常值。子阵为
(2-1-5)
[k rs ] [ Br ]T [ D][B s ]tA
101
二、单元刚度矩阵 1、单元几何矩阵 [ B ] 有了单元的位移模式，利用平面问题的几何方程求得应变分量
0 x x u e e 0 { } [ L][ N ]{} [B ]{} y y v xy y x

三角形常应变单元

第五章P77，求解三角形常应变单元刚度矩阵****************************************************************************** function y=linear_triangle_element_stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)%A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai 0 betaj 0 betam 0;0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammam betam]/(2*A);if p==1D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];else if p==2D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];endendBDDB=D*By=t*A*B'*D*B;******************************************************************************* function y=linear_triangle_assemble(K,k,i,j,m)%K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);K(2*i-1,2*m-1)=K(2*i-1,2*m-1)+k(1,5);K(2*i-1,2*m)=K(2*i-1,2*m)+k(1,6);K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);K(2*i,2*i)=K(2*i,2*i)+k(2,2);K(2*i,2*j-1)=K(2*i,2*j-1)+k(2,3);K(2*i,2*j)=K(2*i,2*j)+k(2,4);K(2*i,2*m-1)=K(2*i,2*m-1)+k(2,5);K(2*i,2*m)=K(2*i,2*m)+k(2,6);K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);K(2*j-1,2*i)=K(2*j-1,2*i)+k(3,2);K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);K(2*j-1,2*m-1)=K(2*j-1,2*m-1)+k(3,5);K(2*j-1,2*m)=K(2*j-1,2*m)+k(3,6);K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);K(2*j,2*i)=K(2*j,2*i)+k(4,2);K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);K(2*j,2*j)=K(2*j,2*j)+k(4,4);K(2*j,2*m-1)=K(2*j,2*m-1)+k(4,5);K(2*j,2*m)=K(2*j,2*m)+k(4,6);K(2*m-1,2*i-1)=K(2*m-1,2*i-1)+k(5,1);K(2*m-1,2*i)=K(2*m-1,2*i)+k(5,2);K(2*m-1,2*j-1)=K(2*m-1,2*j-1)+k(5,3);K(2*m-1,2*j)=K(2*m-1,2*j)+k(5,4);K(2*m-1,2*m-1)=K(2*m-1,2*m-1)+k(5,5);K(2*m-1,2*m)=K(2*m-1,2*m)+k(5,6);K(2*m,2*i-1)=K(2*m,2*i-1)+k(6,1);K(2*m,2*i)=K(2*m,2*i)+k(6,2);K(2*m,2*j-1)=K(2*m,2*j-1)+k(6,3);K(2*m,2*j)=K(2*m,2*j)+k(6,4);K(2*m,2*m-1)=K(2*m,2*m-1)+k(6,5);K(2*m,2*m)=K(2*m,2*m)+k(6,6);y=K;******************************************************************************* function y=linear_tiangle_element_stress(E,NU,xi,yi,xj,yj,xm,ym,p,u)%A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai 0 betaj 0 betam 0;0 gammai 0 gammaj 0 gammam;gammai betai gammaj betaj gammam betam]/(2*A);if p==1D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];else if p==2D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];endendy=D*B*u;******************************************************************************* ******************************************************************************* *************************************************************************clc;E= ;v= ;t= ;k1=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)k2=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)k3=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)k4=linear_triangle_element_stiffness(E,v,t, , , , , , ,1) %(E,v,t,x1,y1,x2,y2, x3,y3,1)K=zeros();%节点坐标个数，节点个数乘以2K=linear_triangle_assemble(K,k1, , , );%节点编号K=linear_triangle_assemble(K,k2, , , );%节点编号K=linear_triangle_assemble(K,k3, , , );%节点编号K=linear_triangle_assemble(K,k4, , , );%节点编号KB=K([ , , , , , , , ,],:) ; %取特定的行5 6 7 8 9 10 11 12k=B(:,[ , , , , , , , ]) %取特定的列5 6 7 8 9 10 11 12f=[ ]’u=k\f %求解位移列向量U=[u(1);0;u(2);u(3);0;0;0;0]%总的位移列向量u1=[U( );U( );U( );U( );U( );U( )];%三角形单元一的位移列向量u2=[U( );U( );U( );U( );U( );U( )]; %三角形单元二的位移列向量u3=[U( );U( );U( );U( );U( );U( )];%三角形单元三的位移列向量u4=[U( );U( );U( );U( );U( );U( )]; %三角形单元四的位移列向量sigma1=linear_tiangle_element_stress(E,v, , , , , , ，1,u1)%单元一的应力大小，(E,v,x1,y1,x2,y2,x3,y3,1,u)sigma1=linear_tiangle_element_stress(E,v, , , , , , ,1,u2) %单元二的应力大小(E,v,x1,y1,x2,y2,x3,y3,1,u)sigma1=linear_tiangle_element_stress(E,v, , , , , , ,1,u3)%单元三的应力大小，(E,v,x1,y1,x2,y2,x3,y3,1,u)sigma1=linear_tiangle_element_stress(E,v, , , , , , ,1,u4) %单元四的应力大小(E,v,x1,y1,x2,y2,x3,y3,1,u)。

三角形常应变单元程序的编制与使用共18页文档

有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。

对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好Array的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。

Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。

本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图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 THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1说明：NNODE—单元结点数，NPION—总结点数， NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向； FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)， DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。

结构力学(5)--结构力学II-1试卷

结构力学II试题1、正误判断题 (每题2分，共16分。

正确的标O，错误的标 )1. 矩阵位移法中各类杆件的单元刚度矩阵均具有对称性和奇异性（）2. 局部坐标系与整体坐标系之间的坐标变换矩阵T是正交矩阵（）3. 简单三角单元的形状函数就是三角单元的面积坐标（）4. 简单三角形单元的位移函数中应该包括常数项（）5. 两类稳定问题的主要区别是：荷载—位移曲线上是否出现分支点（）6. 杆件截面的极限弯矩由极限荷载确定（）7. 由于阻尼的存在，任何振动都不会长期持续下去（）8. 结构计算中，大小、方向随时间变化的荷载必须按动荷载考虑（）二、单项选择题(每题1分，共8分，将选择的字母填在括号内)1. 矩阵位移法的单元刚度方程表示（）A.杆端位移与杆端力之间的关系；B.结点位移与结点力之间的关系；C.结点力与杆端力之间的关系D.结点位移与杆端位移之间的关系2．常应变三角形单元所具有的性质之一是（）A. 各单元之间应力连续；B.单元内任意点各插值函数之和等于1；C. 单元的形状函数完全相等；D.单元的形状函数考虑了面外位移；3. 静力法确定完善体系临界荷载的依据是（）A.未失稳时的静力平衡条件；B. 失稳时的能量守恒条件；C. 失稳时的静力平衡条件；D. 未失稳时的位移协调条件。

4．极限荷载应满足的条件是（）A.单向机构条件 + 内力局限条件；B. 内力局限条件 + 平衡条件；C. 平衡条件 + 单向机构条件；D.平衡条件 + 内力局限条件 + 单向机构条件；5 . 超静定结构的极限荷载的大小（）A. 受截面形状影响；B. 受支座移动影响；C. 受温度变化影响；D. 受材料阻尼影响。

6．某单自由度体系当仅刚度EI扩大4倍时，则周期比原周期（）A、减小1/4；B、减小1/2 ；C、增大1/4；D、增大1/27．忽略直杆的轴向变形，图示结构的动力自由度为（）A 3；B 5；C 4；D 6。

第一节三角形常应变单元

第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。

第一节三角形常应变单元一、结构离散化用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。

(a) (b)图3.1 弹性体和有限元模型将整个结构（平板）划分成有限个三角形。

这样的三角形称为单元（三角形单元）。

三角形单元的顶点取为节点。

3节点三角形单元用边界节点间的直线段来近似板的曲线边界。

这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。

注：1. 全部节点和全部单元一般由1开始按自然顺序编号。

2. 节点编码：总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码每个单元的节点局部码i，j，m和节点总码有一一对应的关系3. 单元间不能有重叠4. 一个单元的任一顶点不许为另一单元任一边的内点5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。

二、位移模式1. 单元节点位移列阵iu图 3.2 平面三角形单元设单元e 的节点号码为i ，j ，m 。

由弹性力学平面可知，单元内任意一点有两个位移分量u ，v ，记为{}Tf u v ⎡⎤⎣⎦=故每个节点也有两个位移分量，因此称节点自由度为2。

3个节点得位移分量分别是 ,,,,,m m i i j j u v u v u v ，用列阵表示为{}i ei i e j j j m m m u v u v u v δδδδ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎩⎭== （3-1）称单元自由度是6。

其中任一子矩阵为{}Ti ii u v δ⎡⎤⎢⎥⎣⎦= （i ，j ，m 轮换）2. 位移模式结构承受载荷以后发生位移，但位移分布事先并不知道。

有限元 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 两个方向的位移,若假定单元位移函数是线性的,则可表示成:(,)123u 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 ed d d d m j j i v u v u v u i {}ii j j m X Y X (2-1-1)Y X Y iej m m F F F F ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪==⎨⎬⎨⎬⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎪⎪⎩⎭确定。

将3个结点坐标(x i,y i ),(x j,y j ),(x m,y m )代入上式得如下两组线性方程:123i i i u x y ααα=++123j j j u x y ααα=++ (a)123m m m u x y ααα=++和546i i i v x y ααα=++546j j j v x y ααα=++ (b)546m m m v x y ααα=++利用线性代数中解方程组的克来姆法则,由(a)可解出待定常数1α 、2α 、3α :11A Aα=22A Aα=33A Aα=式中行列式:1i i i j j j m m m u x y A u x y u x y =2111i i j j m mu y A u y u y =3111i i j jm mx u A x u x u =2111i i j j m mAx y A x y x y ==A 为△ijm 的面积,只要A 不为0,则可由上式解出:11()2m m i ij j a u a u a u A α=++ 21()2m m i ij j bu b u b u A α=++ （C ）31()2m mi i j j c u c u c u A α=++式中:m m i j j a x y x y =- m m j i i a x y x y =- m i j j i a x y x y =-m i j b y y =- m j i b y y =- m i j b y y =- （d ）m i j c x x =- m j i c x x =- m j i c x x =-为了书写方便，可将上式记为: m m i j i a x y x y =- m ij by y =- (,,)i j m m i jc x x =-(,,)i j m 表示按顺序调换下标,即代表采用i,j,m 作轮换的方式便可得到(d)式。

用常应变三角形单元解弹性力学平面问题的程序

用常应变三角形单元解弹性力学平面问题的程序******************************************************************** ANALYSIS PROGTAM OF FINITE ELEMENT METHOD ** FOR PLANE STRESS/STRAIN OF TRIANGULAR ELEMENT ** ----- FEMT3.FOR ----- **------------------------------------------------------------- ** Subroutines: 1-SDATA, 2-STE, 3-ATE, 4-DTE, 5-BTE, 6-STIFF ** 7-EQUPE, 8-INSCD, 9-BGSMT, 10-SIGME ********************************************************************DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)REAL KS(200,100)OPEN(5,FILE='FEMT3.DAT')OPEN(6,FILE='FEMT3.OUT',STATUS='NEW')READ(5,*) NJ,NE,NS,NPJ,IPS(结点、单元、支承、荷载、类型)WRITE(6,*)' FINITE ELEMENT ANALYSIS IN PLANE PROBLEM'WRITE(6,*)' SOURCE DATA OUTPUT'WRITE(6,20) NJ,NE,NS,NPJ,IPS20 FORMAT(4X,'NJ',3X,'NE',3X,'NS',3X,'NPJ',2X,'IPS'/1X,5I5)IF(IPS.EQ.0) WRITE(6,*)' PLANE STRESS PROBLEM'IF(IPS.EQ.1) WRITE(6,*)' PLANE STRAIN PROBLEM'CALL SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,T,V,LND,X,Y,JR,PJ)NJ2=2*NJWRITE(6,50) NJ250 FORMAT(/1X,'DEGREES OF FREEDOM=',I5)WRITE(6,60) NW60 FORMAT(1X,'BAND WIDTH=',I5)CALL STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS)(总刚6)CALL EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P)({P}7)CALL INSCD(NS,NW,NJ2,JR,KS,P)(引入支承条件8)CALL BGSMT(NJ,NJ2,NW,KS,P)（解方程9）CALL SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)（求应力10）CLOSE(5)CLOSE(6)END*--------------------------------------------------------C SUBPROGRAM-1C INPUT STRUCTURAL DATASUBROUTINE SDATA(NJ,NE,NS,NW,NPJ,IPS,E,PR,* T,V,LND,X,Y,JR,PJ)DIMENSION LND(NE,3),X(NJ),Y(NJ),JR(NS,3),PJ(NPJ,3)READ(5,*) E,PR,T,V(弹性模量、泊松比、厚度、容重)WRITE(6,10) E,PR,T,V10 FORMAT(/6X,'E',10X,'PR',9X,'T',9X,'V'/,4F10.2)READ(5,*)((LND(I,J),J=1,3),I=1,NE)（结点编码）WRITE(6,20)20 FORMAT(/1X,'ELEMENT INFORMATION'/3X,'ELEM',3X,* 'I J K'/)WRITE(6,30)(I,(LND(I,J),J=1,3),I=1,NE)30 FORMAT(1X,4I5)READ(5,*)(X(I),Y(I),I=1,NJ)（结点坐标）WRITE(6,40)40 FORMAT(/1X,'COORDINATES OF NODES'/3X,'NODES',* 8X,'X',13X,'Y')WRITE(6,50)(I,X(I),Y(I),I=1,NJ)50 FORMAT(1X,I5,2E15.6)READ(5,*)((JR(I,J),J=1,3),I=1,NS)（约束信息）WRITE(6,60)60 FORMAT(/1X,'CONSTRAINED NODES'/3X,'NODE',3X,'X',4X,'Y') WRITE(6,70)((JR(I,J),J=1,3),I=1,NS)70 FORMAT(1X,3I5)READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息）WRITE(6,80)80 FORMAT(/1X,'LOAD CASES'/3X,'NODE',8X,'X',13X,'Y')WRITE(6,90)((PJ(I,J),J=1,3),I=1,NPJ)90 FORMAT(1X,F5.0,2E15.6)100 NW=0（半带宽）DO 110 IE=1,NEDO 110 I=1,3DO 110 J=1,3IW=IABS(LND(IE,I)-LND(IE,J))IF(NW.LT.IW) THENNW=IWENDIF110 CONTINUENW=(NW+1)*2IF(IPS.NE.0) THENE=E/(1.0-PR*PR)PR=PR/(1.0-PR)ENDIFEND*---------------------------------------------------------C SUBPROGRAM-2C CALCULATE ELEMENT STIFFNESS MATRIXSUBROUTINE STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3)REAL KE(6,6)CALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL DTE(E,PR,D)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,6DO 10 J=1,6KE(I,J)=0.DO 10 K=1,3DO 10 K1=1,310 KE(I,J)=KE(I,J)+B(K,I)*D(K,K1)*B(K1,J)C=AE*TDO 30 I=1,6DO 30 J=1,630 KE(I,J)=KE(I,J)*CEND*------------------------------------------------ C SUBPROGRAM-3C CALCULATE ELEMENT AREASUBROUTINE ATE(IE,NJ,NE,LND,X,Y,AE)DIMENSION LND(NE,3),X(NJ),Y(NJ)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=.5*(XIJ*YIK-XIK*YIJ)END*---------------------------------------------- C SUBPROGRAM-4C CALCULATE ELASTICITY MATRIXSUBROUTINE DTE(E,PR,D)DIMENSION D(3,3)DO 10 I=1,3DO 10 J=1,310 D(I,J)=0.D(1,1)=E/(1.-PR*PR)D(1,2)=E*PR/(1.-PR*PR)D(2,1)=D(1,2)D(2,2)=D(1,1)D(3,3)=.5*E/(1.+PR)END*------------------------------------------------ C SUBPROGRAM-5C CALCULATE MATRIX [B]SUBROUTINE BTE(IE,NJ,NE,LND,X,Y,AE,B)DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6)I=LND(IE,1)J=LND(IE,2)K=LND(IE,3)DO 10 II=1,3DO 10 JJ=1,610 B(II,JJ)=0.B(1,1)=Y(J)-Y(K)B(1,3)=Y(K)-Y(I)B(1,5)=Y(I)-Y(J)B(2,2)=X(K)-X(J)B(2,4)=X(I)-X(K)B(2,6)=X(J)-X(I)B(3,1)=B(2,2)B(3,2)=B(1,1)B(3,3)=B(2,4)B(3,4)=B(1,3)B(3,5)=B(2,6)B(3,6)=B(1,5)DO 20 I1=1,3DO 20 J1=1,620 B(I1,J1)=.5/AE*B(I1,J1)END*------------------------------------------------------- C SUBPROGRAM-6C CALCULATE GLOBAL STIFFNESS MATRIXSUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ)REAL KS(NJ2,NW),KE(6,6)DO 5 I=1,NJ2DO 5 J=1,NW5 KS(I,J)=0.DO 10 IE=1,NECALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)DO 10 I=1,3IZ=LND(IE,I)DO 10 II=1,2IH =2*(I -1)+IIIDH=2*(IZ-1)+IIDO 10 J=1,3JZ=LND(IE,J)DO 10 JJ=1,2L =2*(J -1)+JJIL=2*(JZ-1)+JJIF(IL.GE.IDH) THENIDL=IL-IDH+1KS(IDH,IDL)=KS(IDH,IDL)+KE(IH,L)ENDIF10 CONTINUEEND*-------------------------------------------------------- C SUBPROGRAM-7C CALCULATE NODAL LOAD VECTORSUBROUTINE EQUPE(NJ,NE,NPJ,NJ2,T,V,LND,X,Y,PJ,P) DIMENSION LND(NE,3),X(NJ),Y(NJ),PJ(NPJ,3),P(NJ2) DO 10 I=1,NJ210 P(I)=0.DO 20 I=1,NPJII=PJ(I,1)P(2*II-1)=PJ(I,2)20 P(2*II)=PJ(I,3)30 IF(V.EQ.0.) GOTO 50DO 40 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)PE=-V*AE*T/3.DO 40 I=1,3II=LND(IE,I)40 P(2*II)=P(2*II)+PE50 RETURNEND*---------------------------------------------C SUBPROGRAM-8C INTRODUCE BOUNDARY CONDITIONSUBROUTINE INSCD(NS,NW,NJ2,JR,KS,P)DIMENSION P(NJ2),JR(NS,3)REAL KS(NJ2,NW)DO 30 I=1,NSIR=JR(I,1)DO 30 J=2,3IF(JR(I,J).EQ.0) GOTO 30II=2*IR+J-3KS(II,1)=1.DO 10 JJ=2,NW10 KS(II,JJ)=0.IF(II.GT.NW) JO=NWIF(II.LE.NW) JO=IIDO 20 JJ=2,JO20 KS(II-JJ+1,JJ)=0.P(II)=0.30 CONTINUEEND*-------------------------------------------C SUBPROGRAM-9C SOLVE EQUATIONS BY GS METHODSUBROUTINE BGSMT(NJ,NJ2,NW,KS,P)DIMENSION P(NJ2)REAL KS(NJ2,NW)NJ1=NJ2-1DO 20 K=1,NJ1IF(NJ2.GT.K+NW-1) IM=K+NW-1IF(NJ2.LE.K+NW-1) IM=NJ2K1=K+1DO 20 I=K1,IML=I-K+1C=KS(K,L)/KS(K,1)IW=NW-L+1DO 10 J=1,IWM=J+I-K10 KS(I,J)=KS(I,J)-C*KS(K,M)20 P(I)=P(I)-C*P(K)P(NJ2)=P(NJ2)/KS(NJ2,1)DO 40 I1=1,NJ1I=NJ2-I1IF(NW.GT.NJ2-I+1) JO=NJ2-I+1IF(NW.LE.NJ2-I+1) JO=NWDO 30 J=2,JOK=J+I-130 P(I)=P(I)-KS(I,J)*P(K)40 P(I)=P(I)/KS(I,1)WRITE(6,50)50 FORMAT(/1X,'NODAL DISPLACEMENTS'/3X,* 'NODE',5X,'X-DISP.',8X,'Y-DISP.')DO 60 I=1,NJ60 WRITE(6,70) I,P(2*I-1),P(2*I)70 FORMAT(1X,I5,2E15.6)END*--------------------------------------------------- C SUBPROGRAM-10C CALCULATE ELEMENT STRESS MATRIXSUBROUTINE SIGME(NE,NJ,NJ2,E,PR,LND,X,Y,P)DIMENSION LND(NE,3),X(NJ),Y(NJ),D(3,3),B(3,6), * S(3,6),ST(3),P(NJ2),DE(6)WRITE(6,*)WRITE(6,*)' ELEMENT STRESSES'CALL DTE(E,PR,D)DO 50 IE=1,NECALL ATE(IE,NJ,NE,LND,X,Y,AE)CALL BTE(IE,NJ,NE,LND,X,Y,AE,B)DO 10 I=1,3DO 10 J=1,6S(I,J)=0.DO 10 K=1,310 S(I,J)=S(I,J)+D(I,K)*B(K,J)DO 20 I=1,3DO 20 J=1,2IH=2*(I-1)+JIW=2*(LND(IE,I)-1)+J20 DE(IH)=P(IW)DO 30 I=1,3ST(I)=0.DO 30 J=1,630 ST(I)=ST(I)+S(I,J)*DE(J)SGX=ST(1)SGY=ST(2)TXY=ST(3)ASG=(SGX+SGY)*.5RSG=SQRT(.25*(SGX-SGY)**2+TXY*TXY)SGMA=ASG+RSGSGMI=ASG-RSGIF(SGY.EQ.SGMI) CETA=0.IF(SGY.NE.SGMI) CETA=90.-57.29578*ATAN* (TXY/(SGY-SGMI))50 WRITE(6,60) IE,SGX,SGY,TXY,SGMA,SGMI,CETA60 FORMAT(1X,'ELEMENT NO.=',I4/2X,'SIGX=',E10.4, * 2X,'SIGY=',E10.4,2X,'TXY =',E10.4/2X,'SGMA=', * E10.4,2X,'SGMI=',E10.4,2X,'CETA=',E10.4)END。

8 第四章用常应变三角形单元解力学平面问题 (2)解析

um xm ym
1 um ym
1 xm um
其中
1 xi yi
2 1 x j y j
1 xm ym
(c) (d) (1)
从解析几何可知，式中的就是三角形i、j、m的面积。
为保证求得的面积为正值，节点i、j、m的编排次序必须是逆时针方向，如图1所示。
7. 由单元的节点位移列阵计算单元应力
解出整体结构的节点位移列阵后，再根据单元节点的编号找出对应于单元的位移列阵 e，将 e代入（3-3）式就
可求出各单元的应力分量值。
8. 计算结果输出
求解出整体结构的位移和应力后，可有选择地整理输出某些关键点的位移值和应力值，特别要输出结构的变形图、应力图、应变图、结构仿真变形过程动画图及整体结构的弯矩、剪力图等等。
平面问题可用三角元，四边元等。
例如：
3. 选择单元的位移模式
结构离散化后，要用单元内结点的位移通过插值来获得单元内各点的位移。在有限元法中，通常都是假定单元的位移模式是多项式，一般来说，单元位移多项式的项数应与单元的自由度数相等。它的阶数至少包含常数项和一次项。至于高次项要选取多少项，则应视单元的类型而定。
有限元法的实质是：把有无限个自由度的连续体，理想化为只有有限个自由度的单元集合体，使问题简化为适合于数值解法的结构型问题。
二、经典解与有限元解的区别：
微分经典解法 —— （解析法）
数目增到∞ 大小趋于 0
建立一个描述连续体性质的偏微分方程
有限单元离散化集合
总体分析解
有限元法——连续体——单元——代替原连续体
式中：
Re ke e
（3-4）
——单元刚度矩阵
ke BT DBdxdydz

第一节三角形常应变单元

第一节三角形常应变单元一、结构离散化用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。

(a) (b)将整个结构（平板）划分成有限个三角形。

这样的三角形称为单元（三角形单元）。

三角形单元的顶点取为节点。

3节点三角形单元用边界节点间的直线段来近似板的曲线边界。

这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。

注：1. 全部节点和全部单元一般由1开始按自然顺序编号。

二、位移模式1. 单元节点位移列阵iu设单元e 的节点号码为i ，j ，m 。

由弹性力学平面问题可知，单元内任意一点有两个位移分量u ，v ，记为{}Tf uv ⎡⎤⎣⎦=故每个节点也有两个位移分量，因此称节点自由度为2。

节点i 位移分量记为{}Ti ii u v δ⎡⎤⎣⎦= （i ，j ，m 轮换）则3个节点的位移分量用列阵表示为{}i ei i e j j j m m m u v u v u v δδδδ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎩⎭==（3-1）称为单元节点位移列阵（向量）。

单元自由度是6。

2. 位移模式结构承受载荷以后发生位移，但位移分布事先并不知道。

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

平面三角形单元常应变单元matlaｂ程序的编制
———————————————————————————————— 作者：
———————————————————————————————— 日期：
ﻩ
三角形常应变单元程序的编制与使用
有限元法是求解微分方程边值问题的一种通用数值方法,该方法是一种基于变分法(或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元法中三节点三角形分析结构的步骤如下:
1）整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析,建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效。
POISＳ=fsｃanf（FＰ1，'%f'，1）%泊松比
TＨICＫ=ｆscａnf（ＦＰ1,'%d',1)％厚度
LＮODS=fscanf（ＦＰ1,'%ｄ',［3,NＥLEM］)'%单元定义数组(单元结点号）
功能:运用有限元法中三角形常应变单元解平面问题的计算主程序。
基本思想：单元结点按右手法则顺序编号。
荷载类型：可计算结点荷载。
说明:主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%---－--－-----－---－-----－--－---－－---－--－-----－-－-－－－－-－-－－-－-----------------－----－----－----－----------
glｏｂａｌFＯRCＥFIXＥＤ
globaｌＢMATＸＤMＡTＸSＭATＸAＲＥA

常应变三角形单元

§3 平面问题的有限元分析
§3.1 常应变三角形单元 §3.2 矩形双线性单元 §3.3 有限元分析应注意的问题和结果整理
一.结点的选择和单元划分 1.集中力作用点、分布力突变点、支承点应选作结点。 2.不同厚度、不同材料的部分不应划在同一个单元。 3.应力变化大处单元应密集一些。结点的多少与疏密要考虑计算机的容量和计算精度。 4.单元边界的边长之比应尽可能靠近1。
结点4的应力由结点1、2、3的应力外插得到
F
E
2.两单元平均法三角形单元时，以两相邻单元应力平均值作为边中点的应力近似值。矩形单元时，以两相邻单元公共边两端结点四个应力的平均值作为边中点的应力近似值。对于边界处的结点，同样由内结点结果的外插得到。
总刚需占用的存贮空间为: 16 *14*2=448
三.充分利用结构的对称性 P P P P P
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满意结果，由几何方程求应变，再由物理方程求应力，结果的精度较差。上述三角形单元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有两种:绕结点平均法和两单元平均法。
四.应力结果的整理
位移的计算结果一般比应力、内力结果精度高。位移达到满元为常应力，矩形单元应力线性变化，而工程问题的应力是比较复杂的。为更好地反应实际应力情况，需要对计算结果进行整理。常用处理方法有两种:绕结点平均法和两单元平均法。
1.绕结点平均法以交于同一结点各单元此结点处某应力分量的代数平均值，作为此结点该实际应力的近似值。对于边界处的结点，由内结点结果的外得到。 A B 1 D C 2 3
4
1A B C D E F ( ) 1 1 1 1 1 1 1 6

第6章——常应变三角形单元

收敛——单元尺寸趋于零时，有限元解趋于真解
22:42 有限单元法
崔向阳
12
形函数的性质
当单元的位移函数满足完备性要求时，称单元是完备的（通常较容易满足）。当单元的位移函数满足协调性要求时，称单元是协调的。
当势能泛函中位移函数的导数是2阶时，要求位移函数在单元的交界面上具有C1或更高的连续性，这时构造单元的插值函数往往比较困难。在某些情况下，可以放松对协调性的要求，只要单元能够通过分片试验 (Patch test)，有限元分析的解答仍然可以收敛于正确的解。这种单元称为非协调单元。
(
y1
−
y2
)x
+
(
x2
−
x1
)
y
其中A是三角形的面积
1 A= 1 1
2
x1 x2
y1 y2
1 x3 y3
22:42 有限单元法
崔向阳
5
平面三角形单元
u ( x, y ) = N1u1 + N2u2 + N3u3
同理
v ( x, y ) = N1v1 + N2v2 + N3v3
u(x, y) = N(x, y)de
Ni =1 i
j
i
m Nj =1 j
22:42 有限单元法
崔向阳
i m
j
Nmm =1
9
形函数的性质
在三角形单元的边界ij上任一点（x，y），有:
Ni (x,
y)
=1−
x − xi x j − xi
N j (x,
y)
=
x − xi x j − xi
Nm (x, y) = 0
证

常应变三角形单元

三角形常应变单元程序的编制与使用11页

有限元基础知识归纳

计算力学(有限单元法)第三章重点整理

4平面问题有限元分析

结构力学练习题及答案

9第2章弹性力学平面问题及空间问题有限元

三角形常应变单元

三角形常应变单元程序的编制与使用共18页文档

结构力学(5)--结构力学II-1试卷

第一节三角形常应变单元

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

用常应变三角形单元解弹性力学平面问题的程序

8 第四章 用常应变三角形单元解力学平面问题 (2)解析

第一节三角形常应变单元

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

常应变三角形单元

第6章——常应变三角形单元

8 第四章用常应变三角形单元解力学平面问题 (2)解析