平面三角形单元有限元程序设计.docx
《平面三角形单元有限元程序设计.docx》由会员分享,可在线阅读,更多相关《平面三角形单元有限元程序设计.docx(15页珍藏版)》请在冰豆网上搜索。
![平面三角形单元有限元程序设计.docx](https://file1.bdocx.com/fileroot1/2023-2/22/c0798537-baa0-42fb-8ca0-cf0d6695db45/c0798537-baa0-42fb-8ca0-cf0d6695db451.gif)
平面三角形单元有限元程序设计
、题目如图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)最里层:
区分为第奇/偶数个计算
k1e
66
k1'5656
k2e
66
k2'5656
单元刚度的集成:
kZ
66
kZ'5656
K
ke'ke'
12
kZe'
边界约束的处理:
划0置1法
适用:
这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
做法:
(1)将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1,其他元素改为零;同时将载荷列阵{R}中相应元素用已知位移置换。
◎这样,由该方程求得的此位移值一定等于已知量。
(2)将〔K:
中已知位移相应的列的非主对角成元素也置0,以保持〔K〕的对称性。
◎当然,在已知位移分量不为零的情况下,这样做就改变了方程左端的数值,为保证方程成立,须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。
◎若约束为零位移约束时,此步则可省去。
特点:
1)经以上处理同样可以消除刚性位移(约束足够的前提下),去掉
未知约束反力。
(2)但这种方法不改变方程阶数,利于存贮
(3)不过,
若是要求出约束反力,仍要重新计算各个划去的总刚元
o
程序如下:
变量说明
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
formatshorte
%设定输出类型
clear
%清除内存变量
NELEM=36
%单元个数(单元编码总数)
NPION=28
%结点个数(结点编码总数)
NVFIX=2
%受约束边界点数
NFORCE=1
%结点荷载个数
YOUNG=2e11
%弹性模量
POISS=0.25
%泊松比
%厚度
THICK=0.1
LNODS=[123;245;253;356;
478;485;589;596;
6910;71112;7128;81213;
8139;91314;91410;101415;
111617;111712;121718;121813;
131819;131914;141920;142015;
152021;162223;162317;172324;
172418;182425;182519;251926;
192620;202627;202721;212728]%单元定义数组
(单元结点号)
%相应为单元结点号(编码)、按逆时针顺序输入
COORD=[00;-0.751.5;0.751.5;-1.53;03;
1.53;-2.254.5;-0.754.5;0.754.5;
2.254.5;-36;-1.56;06;1.56;36;
-3.757.5;-2.257.5;-0.757.5;0.757.5;
2.257.5;3.757.5;-4.59;-39;
-1.59;09;1.59;39;4.59]%结点坐标数组
%坐标:
x,y坐标(共NPOIN组)
FORCE=[10-15]%结点力数组(受力结点编号,x方向,y方向)
FIXED=[2211;2811]%约束信息(约束点,x约束,y约束)
%有约束为1,无约
%**********************************************************
%生成单元刚度矩阵并组成总体刚度矩阵
ASTIF=zeros(2*NPION,2*NPION);%生成特定大小总体刚度矩阵并置0
%**********************************************************
fori=1:
NELEM
%生成弹性矩阵D
D=[1
POISS0;
POISS
10;
0
0(1-P0ISS)/2]*Y0UNG心-POISST)
%**********************************************************
%**********************************************************%生成应变矩阵B
forj=0:
2
-COORD(LNODS(i,
b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)(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
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,:
);%临时向量,用来记录当前
单元的节点编号
forj=1:
3
fork=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
%**********************************************************
%将约束信息加入总体刚度矩阵(对角元素改一法)
fori=1:
NVFIX
ifFIXED(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
%**********************************************************
%生成单元刚度矩阵并组成总体刚度矩阵
%**********************************************************ifFIXED(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;%总体荷载向量置
零
fori=1:
NFORCEASLOD((FORCE(i,1)*2-1):
FORCE(i,1)*2)=FORCE(i,2:
3);
end
%**********************************************************
%求解内力
ASDISP=ASTIF\ASLOD'%计算节点位移向量
ELEDISP(1:
6)=0;%当前单元节点位移向量
fori=1:
NELEM
forj=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
Property
halite
Unit
5
铠|ReferencieTemperature
22
6
-诃IsotropicElasbdty
□
7
Derivefrom
Yoirig'smiT
8
Young'sModulus
2E+11
P3工
g
Patsson'sRatio
0.25
10
BiAModulus
r■
1.3333E+11
Pa
11
ShearModulus
fif+iO
Pa
12
Fl泊FieldVariables
-
设置材料参数
U&tflilsotSurtarffBody
¥
GraphicsPrcperti-es
A
-
Definition
Suppressed
NO
StiffnenBehavior
Flexible
CoordmateSystem
MauItCoordinateSystem
ReferenceTemperature
ByEnvironment
Thickness
IGO.mm
ThicknessMt^de
Manual
OffsetTypw
Middle
-
Material
Awignm^nt
StrurturaJStwl
PJortlinearEffects
Yes
V
<
alproject
0或Hodd[A4]
fi-甩Geonejy•--I©SjrfazeBody
二I#去、CoodinabeS巧temw、m、Geb自CosrdirateSystem白®Mest
AJ7nanglesMetiod
日丁白Static5tni(tural(A5}
AnalysisSttirgsE?
画Sokjtkxi(A6)一J]SdutknIrlinratofl
MrrejcciE-eM«W(A4)
三“用LMmeir
L-t]Sjfacexd^
3yCocrdnateSystens
■-了总dobdCoorrtnaie5yst»
T-飾飞时i
:
*助AIT中ng即MEtiod
!
-^EageSiaig
El?
白用bcStndml(A5;,
-(Via^isSeti^j
B?
«jSolitionM
■-_j_Sot^T
k&f"EzzeSizing'-ng
三Scope
乂井ngM^thxl
Scope
SeepingWeticd
Gecrretr>Selectior
Geometty
iBody
Definiitiari
Sdppr«s»d
No
Method
Triangles
|EenientMid^deINodes
UseGlobalSdtnj寸
Z)etaiso:
'AllTrianglesMethod'■Method
Gecmrtry
3Eog«
-Defhition
Supprwsed
hkjrnMrofDivsicns
Nd
Nunb^o2Dvsions
E#c7I0T
-ard
Ed"
NdBias
c_.斤rm
网格划分
'"雪一2曾(wnj
边界约束
添加载荷
A:
Katiettractural
匚qjivdlcri.Sb4
札-qL.ii/^pm卜GlvMirpc]Strns-Top./RntrninUniteMPc
-irrifl:
I
2flT7/ip11;14
CQJ0001S91MMax
QX0C1&93;tjNCUMD
—OJQOCl?
?
^
—DJJOCIMfl?
—099E讥5
—
5JO32C-3
0jC4G7e^
W6^dSMift
A:
MalicStnictunl
Tci柄|口frfoFTTiflran
Type;-Q7fllD©F:
jrQi©*i0i
Unb*:
mm
Time:
'
3017/1/711:
12
Jt6K1e-6Mai
4.1716c&
cxb5D2e-6
3.1237e-6
2XD73e^5
2jOB5Sb-&
k5W4e-6
1.04?
9e-&
5L2145b-7
0Min