八节点等参单元平面有限元

八节点等参单元平面有限元

分析程序

土木工程学院

2011.2

目录

1.概述 (1)

2.编程思想 (2)

2.1.八节点矩形单元介绍 (2)

2.2.有限元分析的模块组织 (5)

3.程序变量及函数说明 (6)

3.1.主要变量说明： (6)

3.2.主要函数说明 (7)

4.程序流程图 (8)

5.程序应用与ANSYS分析的比较 (9)

5.1.问题说明 (9)

5.2.ANSYS分析结果 (10)

5.3.自编程序分析结果 (12)

5.4.结果比较分析 (12)

参考文献 (14)

附录程序源代码 (15)

《计算力学》课程大作业

1.概述

通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。

随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。

C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。

2. 编程思想

本程序采用C 语言编程，编制平面四边形八节点等参元程序，用以求解平面结构问题。程序采用二维等带宽存储整体刚度矩阵，乘大数法引入约束，等带宽高斯消去法求解位移。

在有限元程序中，变量数据需赋值的可分为节点信息，单元信息，载荷信息等。对于一个节点来说，需以下信息：节点编号（整型），节点坐标（实型），节点已知位移（实型），节点载荷（实型），边界条件（实型）和节点温度（实型）等。同样，对于一个单元来说，需以下信息：单元的节点联接信息（整型），单元类型信息（桁架、梁、板、壳等）（整型） ，单元特性信息（厚度、内力矩等）（实型），材料信息（弹性模量，泊松比等）（实型）等。

在FORTRAN 程序中，以上这些变量混合在一起，很难辨认，使程序的可读性不好，如需要进行单元网络的自适应划分，节点及单元的修改将非常困难。在进行C 语言编译过程中，采用结构struct 使每个节点信息存储在一个结构体数组中，提高程序的可读性，使数据结构更趋于合理。

2.1. 八节点矩形单元介绍

八节点矩形单元编号如图 2-1所示

八节点矩形单元的位移函数为：

2

2

2

2

12345678u x y x xy y x y xy

αααααααα=+++++++ (2.1)

2

2

2

2

910111213141516v x y x xy y x y xy αααααααα=+++++++ (2.2)

其形函数为

1122334455667788u N u N u N u N u N u N u N u N u =+++++++ (2.3) 1122334455667788v N v N v N v N v N v N v N v N v =+++++++ (2.4)

式(2.3)和式(2.4)中(,)i i N N εη=，并且采用等参变换，则单元的坐标变换式可取为

1122334455667788

1122334455667788X N x N x N x N x N x N x N x N x Y N y N y N y N y N y N y N y N y =+++++++??

=+++++++?

(2.5)

图 2-1

单元应变矩阵为

{}x y xy u x u y u v y x εεεγ?????

??????????==?????????????

??+??????

(2.6)

式(2.6)一般简写为

{}[]{}B εδ= (2.7) 其中[]B 的子块矩阵为

[]i

i i i i

N x N B y N N y x ??????????

?=????

???

????????

(2.8)

由于i N 是ε、η的函数，在[]i B 中的x 、y 要按照复合函数来求导，即

[]i i

i

i i i N N N x

y x

x J N N N x y y

y εεεη

ηη?????????

????????

?????????????????==?????????????????????????????

????

?

(2.9) 从而有

[]1i i i i N N x

J N N y εη-??????

????

???

???=?????????????

?

?

?

(2.10) 因此，单元应力矩阵为

{}[][]{}D B σδ= (2.11) 单元刚度矩阵为

[][][][]T

e

A

K B D B hdxdy =?? (2.12)

其中积分采用三点高斯积分，

33

11,1

1

1

1

1

(,)()(,)nip

i

j

i j i

i i j i f d d f W

f ξηξηωω

ξηξη--===∑∑∑??

(2.13)

2

nip n =（高斯积分点的总数）

，i ω和j ω或i W 是加权系数，i ξ和j η是单元内的坐标.。对

于三点高斯积分，高斯积分点的位置： 11 5.0ξω==，220.0,8.0ξω==，

33 5.0ξω==。

单元等效节点荷载

{}[]{}T

e

S

P N P ds =? (2.14)

结构刚度矩阵

e

e

K K

=∑ (2.15)

结构结点荷载列阵

e

e

P P

=

∑ (2.16)

注意，对于式(2.15)和式(2.16)中e

∑的理解不是简单的叠加而是集成。 总刚平衡方程

[]{}{}K P δ= (2.17)

从式(2.17)求出

{}[]

{}1

K P δ-= (2.18)

将{}δ回代入式(2.7)和式(2.11)，得到{}ε和{}σ。

2.2.有限元分析的模块组织

一个典型的有限元分析过程主要包括以下几个步骤：

1)读输入数据，定义节点及单元数组。

2)由边界条件计算方程个数，赋值荷载列阵。

3)读入在带状存储的总刚度矩阵中单元和载荷信息。

4)定义总刚度阵数组。

5)组装总刚度阵。

6)解方程组。

其流程图可见图2-2。

图2-2

3.程序变量及函数说明3.1.主要变量说明：

3.2.主要函数说明

4.程序流程图

任何一个有限元程序处理过程，都可以划分为三个阶段：

1)前处理：读入数据和生成数据。

2)处理器：有限元的矩阵计算、组集和求解。

3)后处理：输出节点位移、单元应变等。

为了更清楚地描述程序，我们给出了程序的流程图。

5. 程序应用与ANSYS 分析的比较

为了验证程序的正确性，我们取了一个简单框架结构，分别用自编程序和ANSYS 进行分析和比较。

5.1. 问题说明

图 5-1所示一简支梁，高3m ，长18m ，承受均布荷载210/N m ，10210E Pa =?，

0.167μ=，取1t =m ，作为平面应力问题。

由于结构和荷载对称，只对右边一半进行有限单元计算，如图 5-2所示，而在y 轴各节点布置水平连杆支座。

图 5-1

5.2. ANSYS 分析结果

ANSYS 分析中，采用plane82单元，在板单元上边施加均布荷载210/N m ，在y 轴上的各结点约束UX 方向的自由度，在板单元右下角的结点约束UX 自由度，结点布置、网格划分方案如图 5-3所示。

图 5-3

图 5-4 位移图和图 5-5 各单元X σ图分别给出了结构的位移图和X σ应力云图。

图

5-5 各单元X 图

图 5-4 位移图

从图 5-4 位移图和图 5-5 各单元X σ图可以看到，跨中的位移和正应力最大，与简支梁承受均布荷载，跨中挠度和正应力最大的力学常识相符合，可以初步认为模型是正确的。表格 1给出了0.75x m =的截面上的正应力X σ和表格 2给出了 1.5y m =处各横截面Y 方向位移，其中

X σ的单位为a P ，y δ的单位为m 。

表格 2

5.3. 自编程序分析结果

用自编程序进行分析时，采用了整个结构模型进行计算，其他条件不变，因编程水平有限，在后处理阶段没有给出节点处的位移与应力，而只能得到高斯积分点处的位移与应力，得到积分点处的应力后，根据数值分析相关知识采用了插值进行处理，从而得到0.75x m =的截面上的正应力X σ和 1.5y m =-处各横截面Y 方向位移，分别列出表格如下：

表格 3

表格 4

5.4. 结果分析比较

为了更直观的比较自编程序和ANSYS 的计算结果，将表格 1和表格 3的数据绘入图 5-6，将表格 2和表格 4的数据绘入图 5-7。

图5-6 应力比较

图5-7 位移比较

自编程序所得结果与ANSYS分析结果进行比较发现，应力偏大而位移偏小。分析其中的原因，一方面是编程水平有限，自编程序有很多不完善之处，有很多因素没有考虑（温度、变形、非线性等），所得结果可信度不是很高，只能得到应力以及位移大概的分布情况；另一方面是在后处理阶段，在对高斯积分点结果进行处理时，误差难以避免，还有一方面的原因是在进行求解时保留数据精度不够，造成误差较大，并且进行数值运算时可能会造成大

数吃小数，从而引起结果的差异。

参考文献

[1](美)史密斯(Smith,I.m.)著;王崧等译.有限单元法编程(第三版)[M].北京:电子工业出版

社,2003

[2]王勖成.有限单元法[M].北京:清华大学出版社,2003

[3]宋建辉,涂志刚.MA TLAB语言及其在有限元编程中的应用[J].湛江师范学院学

报.2003.6(24):101-105

[4]郑大玉, 连宏玉, 周跃发. 有限元编程与C语言[J]. 黑龙江商学院学报.1997.3(13):23-28

[5]王伟,刘德富.有限元编程中应用面向对象编程技术的探讨[J].三峡大学学

报.2001.2(23):124-128

[6]汪文立, 吕士俊.二级C语言程序设计教程[M]. 北京:中国水利水电出版社,2006

[7]赵翔龙.FORTRAN 90学习教程[M].北京:北京大学出版社,2002

附录程序源代码

#include

#include

/*The gemotry of the model*/

float mesh[2];

float wide,hight;

float wsize,hsize,young,poiss,thick;

int i,j,k,knelem,npoin,elem[500],ielem;

float coord[2][1],props[3][1],lnods[8][500],ifpre[1]; float posgp[3],weigp[3],estif[136],elcod[2][8];

int npoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb; float gpcod[2][9],bmatx[3][16],dmatx[3][3];

float shape[8],deriv[2][8];

float xjacm[2][2],xjaci[2][2],elcod[2][8];

float cartd[2][8];

float bmatx[3][16],dmatx[3][3];

float nload[1];

float press[3][2],pgash[2],dgash[2];

struct node

{ int nodenu;/*The number of the node*/

float cor[2];/*The coordinate of the node*/

float displacement[2];/*The displacement of the node*/ float load[2]; /*The load of the node*/

float boundary[2];

} node[500];

struct

{ int elementnu;/*The number of element*/

int localnu[8];/*Local number*/

int localcorx[8];

int localcory[8]; /*Local coordinate of node*/

float qx,qy,t;/*The stress and strain*/

}element[500];

void meshing(float wide,float hight,float mesh[2])

{ printf("Please input the meshing size\n");

scanf("%f%f",&wsize,&hsize);

getchar();

mesh[0]=wide/wsize;

mesh[1]=hight/hsize;

}

/*The global coordinate of node*/

void coordinate()

{

int i,j=1;

float x,y;

for(i=1;i<=2*mesh[1]+1;i++)

{x=0;

while(x<=wide)

if(i%2!=0)

{ node[j].cor[1]+=wsize/2;

node[j].cor[2]+=(i-1)*hsize;

j++; }

else

{ node[j].cor[1]+=wsize;

node[j].cor[2]+=(i-1)*hsize;

j++; }

}

}

main()

{ int top[500],bottom[500];/*The number of top and bottom element*/ void input();

void stif();

void jacobi2( );

void load();

void asem();

void solve( );

void stre();

npoin=8;

for(i=1;i<=8;i++)

lnods[i][1]=i;

for(i=1;i<=3;i++)

scanf("%f",&props[i][1]);

printf("The EX is %e\nThe uo is %8f\nThe bt is %8f",props[1][1],props[2][1],props[3][1]); getch();

printf("Please input the gemotry of the model\n");

scanf("%f%f",&wide,&hight);

getchar();

printf("The wide is %0.3f,the hight is %0.3f",wide,hight);

getchar();

meshing(wide,hight,mesh);

printf("The wide size is %f,the hight size is %f\n",mesh[0],mesh[1]);

input();

getchar();

nelem=mesh[0]*mesh[1];

for(i=1;i<=nelem;i++) /*The element number*/

element[i].elementnu=i;

for(i=1;i

npoin+=5;

for(i=1;i

npoin+=3*mesh[0]+2;

printf("%d",npoin);

for(i=1;i<=npoin;i++)

node[i].nodenu=i;

for(i=1;i<=nelem;i++)

for(j=1;j<=8;j++)

element[i].localnu[j]=j;

for (i=1;i<=mesh[0];i++)

bottom[i]=i;

j=1;

for(i=mesh[0]*(mesh[1]-1)+1;i<=nelem;i++)

top[j++]=i;

for(i=1;i

printf("the top number is %d\n",top[i]);

printf("%d\n",element[1].localnu[8]);

coordinate();

stif();

jacobi2( );

load();

asem();

solve( );

stre();

getchar();

}

/*data input*/

void input()

{ int nnode=8;

int notreadchar;

/*input element node number*/

printf("input element node number\n");

for(ielem=1;ielem<=nelem;ielem++)

for(i=1;i<=nnode;i++)

scanf("%d",&lnods[i][ielem]);

/*Input restriction message*/

printf("double-digit is symbol of the restriction,1 representates stable,2 representates gived displacement\n");

i=1;

while(i)

{ scanf("%d",&i);

if(i!=0)

{scanf("%d",&j);

scanf("%d",&ifpre[j]);}

else break;

}

}

/*Element stiffness matrix*/

有限元 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 ?? ?? ???? ???? ??==??????????????????

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

第三章 平面问题的有限元法作业 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 μμμ μμμμμμ---??- - ??????????--?????? -??? ?=????---+ +? ???? ?????????-????+?? ? ? 节点位移

平面问题的有限元法-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 τ

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

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 轴方向。

有限元平面问题

[] ()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 ,,,

平面问题的有限元方法-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.wendangku.net/doc/975382197.html,。祝力学学习顺利！