《有限元基础教程》_【MATLAB算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)

合集下载

如三节点三角形单元

1 Ai 1 x j 2 1 xm
1
x
1 y j ( ai bi x ci y ) 2 ym
y
（7-23）
ai x j ym xm y j bi y j ym ci x j xm
由式（7-23），式（7-21）化为 1 Li (ai bi x ci y ) (i, j , m) 2A
第7章
平面问题高阶单元
7.1 位移模式阶次的选择在前面两章中讨论了平面问题三结点三角形单元，其位移模式的最高阶是坐标 x、 y的一次项。这种位移模式导致单元常应变、常应力特性，单元应变矩阵、应力矩阵、刚度矩阵均为常数矩阵，因此计算非常简单。但这种单元难以反映应力梯度的迅速变化。要想提高计算精度，必须细分网格，增加单元数和点数，因而加大输入数据的工作量提高计算精度的另一条有效途径是采用高阶单元。由于高阶单元的应变、应力不再是常数，因此采用少量单元就
Pmi, Pijd的面积。这三个比值Li, Lj, Lm称为P点的面积坐标。由于则
Ai Ai Am A
Li L j Lm 1
（7-22）
由此可见，P点的三个面积坐标不是独立的。同时，面积
坐标只是用以确定三角形内部某点的位置，因而是一种局部坐标。下面进一步给出面积坐标的几个性质。（2）面积坐标与直角坐标的关系在图7-5中，三角形Pjm的面积为
各待定系数(a1… a8)。将这些系数再代入式（7-1），可得：
u N i ui N j u j N l ul N mum v N i ui N j u j N l ul N mum
式中形函数为：
1 x y N i 1 1 4 a b 1 x y N j 1 1 4 a b 1 x y N l 1 1 4 a b 1 x y N m 1 1 4 a b

《有限元分析》课程作业

《有限元分析》课程作业任课教师：徐亚兰学生姓名：陈新杰学号：班级：1304012时间：2016-01-05一、问题描述及分析问题：如图1所示，有一矩形平板，在右侧受到P=10KN/m 的分布力，材料常数为：弹性模量Pa E 7101⨯=；泊松比3/1=μ；板的厚度为t=；试按平面应力问题利用三角形与矩形单元分别计算各个节点位移及支座反力。

图1 平面矩形结构的有限元分析分析：使用两种方案：一、基于3节点三角形单元的有限元建模，将矩形划分为两个3节点三角形单元；二、基于4节点矩形单元的有限元建模，使用一个4节点矩形单元。

利用MATLAB 软件计算出各要求量，再将两种方案的计算结果进行比较、分析、得出结论。

二、有限元建模及分析1、基于3节点三角形单元的有限元建模及分析（1）结构的离散化与编号如图2所示，将平面矩形结构分为两个3节点三角形单P=10KN/m1m1m元。

单元①三个节点的编号为1，2，4，单元②三个节点的编号为3，4，2，各个节点的位置坐标为(),,1,2,3,4i i x y i =，各个节点的位移（分别沿x 方向和y 方向）为(),,1,2,3,4i i u v i =。

图2 方案一：使用两个3节点三角形单元（2）各单元的刚度矩阵及刚度方程 a.单元的几何和节点描述单元①有6个节点位移自由度（DOF ）。

将所有节点上的位移组成一个列阵，记作(1)q ；同样，将所有节点上的各个力也组成一个列阵，记作(1)F ，则有(1)112244,,,,,)q u v u v u v =（(1)112244(,,,,,)x y x y x y F F F F F F F =同理，对于单元②，有(2)334422,,,,,)q u v u v u v =（1234X y ①②(2)334422(,,,,,)x y x y x y F F F F F F F =b.单元的位移场描述对于单元①，设位移函数012012(,)(,)u x y a a x a y v x y b b x b y ⎫=++⎪⎬=++⎪⎭(1-1)由节点条件，在,i i x x y y ==处，有(,)(,)i i i i i i u x y u v x y v =⎫⎬=⎭1,2,4i = (1-2) 将式(1-1)代入节点条件式(1-2)中，可求出式（1-1）中待定系数，即011122211223444411()22u x y a u x y a u a u a u AAu x y ==++ (1-3) 11122112234441111()221u y a u y b u b u b u AAu y ==++ (1-4) 21122112234441111()221x u a x u c u c u c u AAx u ==++ (1-5) 01122341()2b a v a v a v A =++(1-6) 11122341()2b b v b v b v A =++(1-7) 21122341()2b c v c v c v A =++(1-8)在式(1-3)~式(1-8)中1122123441111()221x y A x y a a a x y ==++ (1-9)2212442442124421244(1,2,3)1111x y a x y x y x y y b y y y x c x x x ⎫==-⎪⎪⎪⎪=-=-⎬⎪⎪⎪==-+⎪⎭ (1-10) 上式中的符号(1,2,3)表示下标轮换，如12,23,31→→→同时更换。

有限元程序设计--第五章线性三角形单元

03:25
16
单元应变和应力矩阵
由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常数矩阵。也就是说，三角形三节点单元内的应力分量也是常量。当然，相邻单元的E，， A和bi、ci(i，j， m)一般不完全相同，因而具有不同的应力，这就造成在相邻单元的公共边上存在着应力突变现象。但是随着网格的细分，这种突变将会迅速减小。
b2 B2 0 c2
b1 , c1 , b2 , c2 , b3 , c3 与x、y无关，都是常量，因此B矩
阵也是常量。单元中任一点的应变分量是 B矩阵与单元节点位
移的乘积，因而也都是常量。因此，这种单元被称为常应变单
元。
03:25
15
单元应变和应力矩阵
平面应力： x σ ( x, y ) y
03:25
12
形函数的性质
当单元的位移函数满足完备性要求时，称单元是完备的（通
常较容易满足）。当单元的位移函数满足协调性要求时，称单元是协调的。
当势能泛函中位移函数的导数是2阶时，要求位移函数在单元
的交界面上具有C1或更高的连续性，这时构造单元的插值函数往往比较困难。在某些情况下，可以放松对协调性的要求，只要单元能够通过分片试验 (Patch test)，有限元分析的解答仍然可以收敛于正确的解。这种单元称为非协调单元。
3
平面三角形单元
显然，三角形三个节点的的位移可由下列方程给出，在各节点上的水平位移方程为： u1=1+2 x1 +3 y1 u2=1+2 x2 +3 y2 u3=1+2 x3 +3 y3
1 1 x1 解出 2 1 x2 1 x 3 3 y1 y2 y3

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

基于Matlab说话的按平面三角形单元划分的构造有限元程序设计专业：建筑与土木匠程班级：建工研12-2姓名：韩志强学号：471220580基于Matlab说话的按平面三角形单元划分构造有限元程序设计一、有限单元发及Matlab说话概述1. 有限单元法跟着现代工业.临盆技巧的成长,不竭请求设计高质量.高程度的大型.庞杂和周详的机械及工程构造.为此目标,人们必须预先经由过程有用的盘算手腕,确实的猜测即将诞生的机械和工程构造,在将来工作时所产生的应力.应变和位移是以,须要追求一种简略而又准确的数值剖析办法.有限单元法恰是顺应这种请求而产生和成长起来的一种十分有用的数值盘算办法.有限元法把一个庞杂的构造分化成相对简略的“单元”,各单元之间经由过程结点互相衔接.单元内的物理量由单元结点上的物理量按必定的假设内插得到,如许就把一个庞杂构造从无穷多个自由度简化为有限个单元构成的构造.我们只要剖析每个单元的力学特征,然后按照有限元法的规矩把这些单元“拼装”成整体,就可以或许得到整体构造的力学特征.有限单元法根本步调如下：(1)构造离散：构造离散就是树立构造的有限元模子,又称为网格划分或单元划分,即将构造离散为由有限个单元构成的有限元模子.在该步调中,须要依据构造的几何特征.载荷情形等肯定单元体内随意率性一点的位移插值函数.(2)单元剖析：依据弹性力学的几何方程以及物理方程肯定单元的刚度矩阵.(3)整体剖析:把各个单元按本来的构造从新衔接起来,并在单元刚度矩阵的基本上肯定构造的总刚度矩阵,形成如下式所示的整体有限元线性方程：式中阵.(4)载荷移置:依据静力等效道理,将载荷移置到响应的节点上,形成节点载荷矩阵.(5)鸿沟前提处理:对式①所示的有限元线性方程进行鸿沟前提处理.(6)求解线性方程:求解式①所示的有限元线性方程,得到节点的位移.在该步调中,如有限元模子的节点越多,则线性方程的数目就越多,随之有限元剖析的盘算量也将越大.(7)求解单元应力及应变依据求出的节点位移求解单元的应力和应变.(8)成果处理与显示进入有限元剖析的后处理部分,对盘算出来的成果进行加工处理,并以各类情势将盘算成果显示出.在用有限元法进行构造剖析时,将会碰到大量的数值盘算,因而在实用上是必定要借助于盘算机和有限元程序,才干完成这些庞杂而沉重的数值盘算工作.而Matlab是当今国际科学界最具影响力和活气的软件.它来源于矩阵运算,并已经成长成一种高度集成的盘算机说话.它供给了壮大的科学盘算,灵巧的程序设计流程,高质量的图形可视化与界面设计,便捷的与其他程序和说话接口的功效.Matlab在列国高校与研讨单位起侧重大的感化.“工欲善其事,必先利其器”.假如有一种十分有用的对象能解决在教授教养与研讨中碰到的问题,那么Matlab说话恰是如许的一种对象.它可以将运用者从繁琐.无谓的底层编程中解放出来,把有限的珍贵时光更多地花在解决问题中,如许无疑会进步工作效力. 今朝,Matlab已经成为国际上最风行的科学与工程盘算的软件对象,如今的Matlab已经不但仅是一个“矩阵试验室”了,它已经成为了一种具有普遍运用远景的全新的盘算机高等编程说话了,有人称它为“第四代”盘算机说话,它在国表里高校和研讨部分正扮演侧重要的脚色.Matlab说话的功效也越来越壮大,不竭顺应新的请求提出新的解决办法.可以预感,在科学运算.主动掌握与科学画图范畴Matlab说话将长期保持其举世无双的地位.为此,本例采取Matlab说话编程,以运用其快捷壮大的矩阵数值盘算功效.二、问题描写一矩形薄板,一边固定,推却150kN分散荷载,构造简图及按平面三角形单元划分的有限元模子图如下所示.薄板厚在本例中,所取构造模子及数据重要用于程序设计理论剖析,与工程现实无关.三、参数输入：单元个数NELEM=4节点个数NNODE=6受束缚鸿沟点数NVFIX=2节点荷载个数NFORCE=1弹性模量YOUNG=2e8泊松比POISS=0.2厚度THICK=0.002单元节点编码数组节点坐标数组节点力数组FORCE =[4 0 -150]束缚信息数组以上数值数据为程序运行前输入的初始数据,存为“”文本格局,同时必须放在Matlab工作目次下,路径不合错误程序不克不及主动读取指定初始文件,运行出错.初始数据文本文件输入格局如下图：四、Matlab说话程序源代码：1.程序中变量解释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 数据文件指针2.Matlab说话程序代码%****************************************************************************** %初始化及数据挪用format short e %设定输出类型clear %消除内存变量FP1=fopen('','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])' %结点力数组（受力结点编号,x偏向,y偏向）FIXED=fscanf(FP1,'%d',[3,NVFIX])' %束缚信息（束缚点,x束缚,y束缚）%有束缚为1,无束缚为0%***************************************************************************** %生成单元刚度矩阵并构成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0%***************************************************************************** for i=1:NELEM%生成弹性矩阵DD= [1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)%***************************************************************************** %盘算当前单元的面积AA=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%***************************************************************************** %生成应变矩阵Bfor j=0:2b(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);endB=[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*BS=D*B;ESTIF=B'*S*THICK*A; %求解单元刚度矩阵a=LNODS(i,:); %暂时向量,用来记载当前单元的节点编号for j=1:3for k=1:3ASTIF((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); %依据节点编号对应关系将单元刚度分块叠加到总刚%度矩阵中endendend%***************************************************************************** %将束缚信息参加总体刚度矩阵（对角元素改一法）for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(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;%对角元素为1endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0; %一列为零ASTIF(FIXED(i,1)*2,:)=0; %一行动零ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1endend%***************************************************************************** %生成荷载向量ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%***************************************************************************** %求解内力ASDISP=ASTIF\ASLOD' %盘算节点位移向量ELEDISP(1:6)=0;%当前单元节点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%掏出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP' %求内力endfclose(FP1); %封闭数据文件五、程序运行成果NELEM =4NPION =6NVFIX =2NFORCE =1YOUNG =200000000POISS =THICK =LNODS =1 2 62 3 42 4 52 5 6COORD =0 01 02 02 11 10 1FORCE =4 0 -150FIXED =1 1 16 1 1D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 00 0 8.3333e+007A =D =2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007A =ASDISP =（解释：以上12个ASDISP输出成果数据暗示节点位移向量,即各节点分离在X,Y偏向上产生的位移.）i =1STRESS =-1.6793e+005-3.3586e+004-1.3207e+005i =2STRESS =-5.5503e+004-5.5503e+004-5.5503e+004i =3STRESS =5.5503e+004-4.9526e+004-9.4497e+004i =4STRESS =1.6793e+005-2.7040e+004-1.7932e+004（解释：以上四组STRESS输出成果数据暗示单元应力SX,SY,SXY,i为单元号.）六、ANSYS建模比较剖析运用ANSYS有限元剖析软件,完整按照前面Matlab程序设计的各项前提参数以及单元划分方法树立ANSYS模子,其求解成果与以上程序盘算成果比较,验证程序是否准确.ANSYS模子节点单元如下图所示：ANSYS模子求解变形图如下所示：ANSYS求解节点位移成果列表显示如下：（单位：m）PRINT U NODAL SOLUTION PER NODE***** POST1 NODAL DEGREE OF FREEDOM LISTING *****THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATESYSTEMNODE UX UY UZ USUM1 0.0000 0.0000 0.0000 0.00006 0.0000 0.0000 0.0000 0.0000MAXIMUM ABSOLUTE VALUESNODE 4 4 0 4ANSYS求解单元应力成果列表显示如下：（单位：KN/m2）PRINT S ELEMENT SOLUTION PER ELEMENT***** POST1 ELEMENT NODAL STRESS LISTING *****THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATESELEMENT= 1 PLANE182NODE SX SY SZ SXY SYZ1 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.00002 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182NODE SX SY SZ SXY SYZ2 -55503. -55503. 0.0000 -55503. 0.00003 -55503. -55503. 0.0000 -55503. 0.00004 -55503. -55503. 0.0000 -55503. 0.0000 ELEMENT= 3 PLANE182NODE SX SY SZ SXY SYZ2 55503. -49526. 0.0000 -94497. 0.00004 55503. -49526. 0.0000 -94497. 0.00005 55503. -49526. 0.0000 -94497. 0.0000 ELEMENT= 4 PLANE182NODE SX SY SZ SXY SYZ2 0.16793E+06 -27040. 0.0000 -17932. 0.00005 0.16793E+06 -27040. 0.0000 -17932. 0.00006 0.16793E+06 -27040. 0.0000 -17932. 0.0000结论经由过程比较Matlab说话设计程序运行成果和ANSYS建模剖析成果可知,两种方法算出的成果完整一致,说程序设计准确.所以本程序实用于按三角形单元划分的平面构造有限元剖析.format short eclearFP1=fopen('LinearTriangleElement of Thin plate.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])FORCE=fscanf(FP1,'%f',[3,NFORCE])FIXED=fscanf(FP1,'%d',[3,NVFIX])ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS^2)*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]%盘算当前单元的面积AA=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%生成应变矩阵Bfor j=0:2b(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); endB=[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*BS=D*B;ESTIF=B'*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF((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);endendend%将束缚信息参加总体刚度矩阵for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(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;endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIF\ASLOD'ELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);endiSTRESS=D*B1(:, :, i)*ELEDISP'endfclose(FP1);format short eclearFP1=fopen('LinearTriangleElement of Thin plate.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',[NELEM,3])COORD=fscanf(FP1,'%f',[NPION,2])FORCE=fscanf(FP1,'%f',[NFORCE,3])FIXED=fscanf(FP1,'%d',[NVFIX,3])ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 for i=1:NELEM%生成弹性矩阵DD= YOUNG/(1-POISS^2)*[1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]%盘算当前单元的面积AA=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%生成应变矩阵Bfor j=0:2b(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); endB=[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*BS=D*B;ESTIF=B'*S*THICK*A;a=LNODS(i,:);for j=1:3for k=1:3ASTIF((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);endendend%将束缚信息参加总体刚度矩阵for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(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;endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0;ASTIF(FIXED(i,1)*2,:)=0;ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1;endend%生成荷载向量ASLOD(1:2*NPION)=0;for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3)end%求解内力ASDISP=ASTIF\ASLOD'ELEDISP(1:6)=0;for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); endiSTRESS=D*B1(:, :, i)*ELEDISP'endfclose(FP1);。

基于MTLAB的3节点三角形单元求解平面弹性问题

本文部分内容来自网络整理，本司不为其真实性负责，如有异议或侵权请及时联系，本司将立即删除！== 本文为word格式，下载后可方便编辑和修改！ ==深知兵真爱兵,打牢官兵思想篇一：对战士的爱更深沉些对战士的爱更深沉些、对战士关心更长远些一直以来，总队医院党委始终把促进战士成长、成才作为提高部队战斗力的重要举措。

针对战士整体素质高但个体差异较大、所学知识多但弄懂搞精的少的实际，医院党委想方设法为战士学习成才创造条件，尽力满足战士求知成才愿望。

使医院开展的“深知兵、真爱兵、带好兵”活动迸发出更大活力。

一、在“深知”上下功夫。

知兵是爱兵的基础。

总队医院党委不仅要求机关干部和带兵干部在了解大量的基本数据，掌握思想变化的“活情况”的同时，而且深层次了解战士能干什么，给每名战士“量体裁衣”，为每个战士建立成才档案，每年在新战士下连后，总队医院都会开展“走近士兵、研究士兵”活动。

主动摸清战士需求，要对战士的长远发展负责，引导他们把个人追求与部队建设有机统一起来，在岗位上成才，在成才中进步。

为战士合理制定成才计划，使战士的合理愿望切合实际。

设身处地的关爱，拉近了官兵心与心之间的距离，战士有话愿向干部讲，有事愿靠组织办，有苦愿向骨干诉，有乐愿与大家享的氛围已经形成，上下和谐，官兵关系顺畅。

二、在“真爱”上做文章。

总队医院积极为战士成才创造条件，医院党委投资100多万元为中队建立电脑网络学习室，与驻地4家大型企业建立了战士成才协作机制。

总队医院充分运用医院自身人才技术队伍力量，对战士进行计算机操作、医学理疗、超声波技术、放射技术等专项技能培训。

为了使培训取得实效，使战士的计算机技能得到社会认可，医院还主动邀请、积极协调地方人力资源部门专业人员到部队授课，并对参加培训的战士进行了职业技能鉴定考核。

目前，参加各类培训的战士陆续如愿获得各种职业资格证书。

同时，政治处鼓励战士参加军地自考、函授学习，并借助各种优势，邀请地方院校教授来部队作专题知识讲座，为战士学习成才搭建平台，调动了战士学习成才积极性。

matlab有限元三单元编程

matlab有限元三单元编程MATLAB是一种功能强大的编程语言和环境，广泛用于工程和科学领域的数值计算、数据分析和可视化。

在有限元分析中，MATLAB的强大功能和直观的语法使其成为一个理想的选择。

在本文中，我们将讨论MATLAB中有限元三单元的编程方法和实践。

有限元分析是一种数值方法，用于解决连续介质的力学问题。

它将一个复杂的结构分解成更简单的有限元单元，然后通过求解线性代数方程组来得到结构的应力和位移解。

在有限元分析中，三角形和四边形单元是最常用的有限元单元之一。

本文将重点讨论三角形单元的编程实现。

首先，我们需要定义一个三角形单元的几何信息。

在三角形单元中，我们有三个顶点，每个顶点有两个坐标。

我们可以使用一个3x2的矩阵来表示这些坐标。

例如，定义三角形ABC的顶点坐标矩阵为P：P = [x_A, y_A;x_B, y_B;x_C, y_C];接下来，我们需要定义三角形单元的连接性信息。

在MATLAB中，我们可以使用一个3x1的向量来表示三个顶点的连接性。

例如，定义三角形ABC的连接性向量为E：E = [1;2;3];然后，我们可以定义三角形单元材料属性和载荷信息。

这些信息包括杨氏模量、泊松比和外力向量。

我们可以将这些信息存储在一个结构体中，例如：properties.E = 210e9; % 杨氏模量properties.v = 0.3; % 泊松比properties.f = [0; -1000; 0]; % 外力向量接下来，我们可以使用这些信息来计算三角形单元的刚度矩阵和力向量。

刚度矩阵是一个3x3的矩阵，力向量是一个3x1的向量。

我们可以使用以下公式来计算它们：K = (E/(1-v^2)) * [1 v 0;v 1 0;0 0 (1-v)/2] * A;f = (A/3) * [1; 1; 1] * properties.f;其中，A是三角形的面积。

我们可以使用以下公式来计算它：A = abs(det([1, 1, 1;x_A, x_B, x_C;y_A, y_B, y_C])) / 2;最后，我们可以使用这些信息来求解三角形单元的位移解。

有限元分析基础教程

有限元分析基础教程前言有限元分析已经在教学、科研以及工程应用中成为重要而又普及的数值分析方法和工具；该基础教程力求提供具备现代特色的实用教程。

在教材的内容体系上综合考虑有限元方法的力学分析原理、建模技巧、应用领域、软件平台、实例分析这几个方面，按照教科书的方式深入浅出地叙述有限元方法，并体现出有限元原理“在使用中学习，在学习中使用”的交互式特点，在介绍每一种单元的同时，提供完整的典型推导实例、MATLAB实际编程以及ANSYS应用数值算例，并且给出的各种类型的算例都具有较好的前后对应性，使学员在学习分析原理的同时，也进行实际编程和有限元分析软件的操作，经历实例建模、求解、分析和结果评判的全过程，在实践的基础上深刻理解和掌握有限元分析方法。

一本基础教材应该在培养学员掌握坚实的基础理论、系统的专业知识方面发挥作用，因此，教材不但要提供系统的、具有一定深度的基础理论，还要介绍相关的应用领域，以给学员进一步学习提供扩展空间，本教程正是按照这一思路进行设计的；全书的内容包括两个部分，共分9章；第一部分为有限元分析基本原理，包括第1章至第5章，内容有：绪论、有限元分析过程的概要、杆梁结构分析的有限元方法、连续体结构分析的有限元方法、有限元分析中的若干问题讨论；第二部分为有限元分析的典型应用领域，包括第6章至第9章，内容有：静力结构的有限元分析、结构振动的有限元分析、传热过程的有限元分析、弹塑性材料的有限元分析。

在基本原理方面，以基本变量、基本方程、求解原理、单元构建等一系列规范的方式进行介绍；在阐述有限元分析与应用方面，采用典型例题、MATLAB程序及算例、ANSYS算例的方式，以体现出分析建模的不同阶段和层次，引导学员领会有限元方法的实质，还提供有大量的练习题。

本教程的重点是强调有限元方法的实质理解和融会贯通，力求精而透，强调学员综合能力(掌握和应用有限元方法)的培养，为学员亲自参与建模、以及使用先进的有限元软件平台提供较好的素材；同时，给学员进一步学习提供新的空间。

matlab有限元三角形单元编程

matlab有限元三角形单元编程
在MATLAB中进行有限元分析，可以使用其自带的有限元分析工具箱（FEATool）进行编程。

以下是一个简单的例子，演示如何使用三角形单元进行有限元分析：
1. 打开MATLAB，进入FEATool环境。

2. 创建新的有限元模型。

选择“Model”菜单下的“New Model”选项，进入“Model Builder”界面。

3. 在“Model Builder”界面中，选择“2D Triangle”单元类型，并在绘图区域中绘制出三角形网格。

4. 在“Model Builder”界面中，设置材料属性、边界条件和载荷等参数。

5. 运行有限元分析。

选择“Model”菜单下的“Solve”选项，进行有限元求解。

6. 查看结果。

选择“Model”菜单下的“Results”选项，可以查看位移、应力、应变等结果。

以上是一个简单的例子，演示了如何使用三角形单元进行有限元分析。

在实际应用中，还需要根据具体问题进行详细的建模和计算。

matlab有限元三角形单元程序

matlab有限元三角形单元程序以下是一个简单的 MATLAB 有限元三角形单元程序的示例：```matlab% 定义模型参数E = 1000; % 弹性模量nu = 0.3; % 泊松比thickness = 1; % 板的厚度% 定义节点坐标nodes = [0, 0; 1, 0; 1, 1; 0, 1; 0.5, 0.5];% 定义单元连接关系connectivity = [1, 2, 5; 2, 3, 5; 3, 4, 5; 4, 1, 5];% 计算总节点数和总单元数numNodes = size(nodes, 1);numElements = size(connectivity, 1);% 初始化全局刚度矩阵和载荷向量K = zeros(numNodes);F = zeros(numNodes, 1);% 循环遍历每个单元for i = 1:numElements% 查找当前单元的节点编号nodesIndex = connectivity(i, :);% 根据节点编号从全局坐标矩阵中取出节点坐标coordinates = nodes(nodesIndex, :);% 计算当前单元的局部刚度矩阵localK = calculateLocalStiffness(E, nu, thickness, coordinates);% 组装局部刚度矩阵到全局刚度矩阵中K(nodesIndex, nodesIndex) = K(nodesIndex, nodesIndex) + localK;% 计算当前单元的局部载荷向量localF = calculateLocalLoad(thickness, coordinates);% 组装局部载荷向量到全局载荷向量中F(nodesIndex) = F(nodesIndex) + localF;end% 边界条件：节点1固定K(1, :) = 0;K(1, 1) = 1;F(1) = 0;% 解线性方程组U = K \ F;% 输出位移结果disp('节点位移：');disp(U);% 计算应力结果stress = calculateStress(E, nu, thickness, nodes, connectivity, U);% 输出应力结果disp('节点应力：');disp(stress);% 计算局部刚度矩阵的函数function localK = calculateLocalStiffness(E, nu, thickness, coordinates)% 计算单元的雅可比矩阵J = (1/2) * [coordinates(2,1)-coordinates(1,1), coordinates(3,1)-coordinates(1,1);coordinates(2,2)-coordinates(1,2), coordinates(3,2)-coordinates(1,2)];% 计算雅可比矩阵的逆矩阵invJ = inv(J);% 计算单元刚度矩阵B = invJ * [-1, 1, 0; -1, 0, 1];D = (E/(1-nu^2)) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; localK = thickness * abs(det(J)) * (B' * D * B);end% 计算局部载荷向量的函数function localF = calculateLocalLoad(thickness, coordinates) localF = zeros(3, 1);midPoint = [sum(coordinates(:,1))/3,sum(coordinates(:,2))/3];localF(3) = thickness * 1 * det([coordinates(1,:); coordinates(2,:); midPoint]);end% 计算各节点应力的函数function stress = calculateStress(E, nu, thickness, nodes, connectivity, U)stress = zeros(size(nodes, 1), 3);for i = 1:size(connectivity, 1)nodesIndex = connectivity(i, :);coordinates = nodes(nodesIndex, :);Ke = calculateLocalStiffness(E, nu, thickness, coordinates);Ue = U(nodesIndex);stress(nodesIndex, :) = stress(nodesIndex, :) + (Ke * Ue)';endstress = stress / thickness;end```这个程序实现了一个简单的平面三角形单元有限元分析，包括定义节点坐标和单元连接关系、计算全局刚度矩阵和载荷向量、施加边界条件、解线性方程组、计算节点位移和应力等。

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

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

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

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

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

本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。

有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。

2）单元分析，建立单元刚度矩阵。

3）整体分析，建立总刚矩阵。

4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。

6）解方程，求出节点位移。

7）求出各单元的单元应力。

8）计算结果整理。

计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。

图1 程序流程图1.1 程序说明%******************************************************************* % 三角形常应变单元求解结构主程序%******************************************************************* ●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。

《有限元基础教程》_【MATLAB算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)

基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)如图4-20所示为一矩形薄平板，在右端部受集中力100 000F N =作用，材料常数为：弹性模量7110E Pa =⨯，泊松比13μ=，板的厚度0.1t m =。

基于MA TLAB 平台求解该结构的节点位移、支反力以及单元应力。

图4-20解答：对该问题进行有限元分析的过程如下。

（1）结构的离散化与编号将结构离散为二个3节点三角形单元，单元编号及节点编号如图4-20(b)所示。

（2）计算各单元的刚度矩阵(以国际标准单位)首先在MA TLAB 环境下，输入弹性模量E 、泊松比NU 、薄板厚度t 和平面应力问题性质指示参数ID ，然后针对单元1和单元2，分别两次调用函数Triangle2D3Node_Stiffness ，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)。

>> E=1e7;>> NU=1/3;>> t=0.1;>> ID=1;>> k1=Triangle2D3Node_Stiffness(E,NU,t,2,0,0,1,0,0,ID)k1 = 1.0e+006 *0.2813 0 0 0.1875 -0.2813 -0.18750 0.0938 0.1875 0 -0.1875 -0.09380 0.1875 0.3750 0 -0.3750 -0.18750.1875 0 0 1.1250 -0.1875 -1.1250-0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750-0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188>>k2=Triangle2D3Node_Stiffness(E,NU,t,0,1,2,0,2,1,ID)k2 = 1.0e+006 *0.2813 0 0 0.1875 -0.2813 -0.18750 0.0938 0.1875 0 -0.1875 -0.09380 0.1875 0.3750 0 -0.3750 -0.18750.1875 0 0 1.1250 -0.1875 -1.1250-0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750-0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188（3）建立整体刚度方程由于该结构共有4个节点，则总共的自由度数为8，因此，结构总的刚度矩阵为KK (8×8)，先对KK 清零，然后两次调用函数Triangle2D3Node_Assembly 进行刚度矩阵的组装。

3结点三角形单元有限元程序MATLAB语言

3结点三角形单元有限元程序（MATLAB语言）该程序包括以下6个部分：1．主程序tri_fem：用于数据的录入和其他程序的调用；2．总刚程序Kf：计算结构的总体刚度；3．各结点位移求解程序xf：求解各结点的位移；4．线性方程组求解程序Jordan：Gauss-Jordan法求解非约束结点的位移；5．应力应变程序ss：由各结点位移求解各单元内的三个结点的应力stress和应变strain；6．数据录入程序input：录入材料、几何尺寸、单元编号和结点编号、位移约束和已知载荷等。

以课本P25页例2.2为例，其input程序为function [E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input()E=2.1e11; v=1/3; t=1; %杨氏模量Pa，泊松比，厚度EN=2; %单元数ecode=[3 1 2; %单元编号单元1 3-1-2;单元2 1-3-41 3 4];NN=4; %结点数node=[0 0; %各结点坐标2 0;2 1;0 1];RN=2; %被约束的位移数RC=[1 4]; %被约束的结点PN=2; %有载荷的结点数%PC（1）表示有载荷的结点，PC(2)表示各结点的力，PC（3）表示载荷方向，0为x方向，1为y方向PC=[2 3;-1/2 -1/2;1 1]; %结点2、3分别有y负方向上-1/2N的力作用在matlab环境下，输入则程序运行结果如下：该程序求解的结点位移结果和结点应力结果与课本给出的结果一致。

附录：%%-------------平面三角形单元有限元法---------------------%%function [x strain stress]=tri_fem()[E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input; %调入已定材料、几何尺寸以及单元和结点编号及约束和载荷分布[n m]=size(ecode);if EN~=nerror('Wrong elementnumber EN or wrong elementcode ecode!');return;end[n m]=size(node);if NN~=nerror('Wrong nodenumber NN or wrong node-coordinate node!');return;ende=zeros(EN,6);A=zeros(EN,1); %面积for i=1:ENe(i,:)=[node(ecode(i,1),:),node(ecode(i,2),:),node(ecode(i,3),:)]; %各三角形单元的节点坐标D=[1,node(ecode(i,1),:)1,node(ecode(i,2),:)1,node(ecode(i,3),:)];A(i,1)=1/2*det(D);end%% 形成单元参数b=zeros(EN,3);c=zeros(EN,3); %各单元参数初始化for i=1:ENb(i,1)=e(i,4)-e(i,6); b(i,2)=e(i,6)-e(i,2); b(i,3)=e(i,2)-e(i,4);c(i,1)=-e(i,3)+e(i,5); c(i,2)=-e(i,5)+e(i,1); c(i,3)=-e(i,1)+e(i,3);end%% 求得总刚，并引入约束和载荷求得各结点位移K=Kf(E,v,t,EN,ecode,NN,A,b,c); %调用函数Kf，求得结构的总体刚度矩阵x=xf(NN,RN,RC,PN,PC,K); %调用函数xf,求得各结点位移%% 求解应力应变[strain stress]=ss(E,v,EN,ecode,A,b,c,x);%% 单元刚度矩阵与结构刚度矩阵function K=Kf(E,v,t,EN,ecode,NN,A,b,c)Ke=zeros(6,6); %单元的刚度矩阵，初始为6*6阶零矩阵K=zeros(NN*2,NN*2); %结构的总体刚度矩阵，初始为零矩阵for m=1:1:EN %m为单元号for i=1:1:3for j=1:1:3Ke(2*i-1,2*j-1)=b(m,i)*b(m,j)+(1-v)*c(m,i)*c(m,j)/2;Ke(2*i-1,2*j)=v*c(m,i)*b(m,j)+(1-v)*b(m,i)*c(m,j)/2;Ke(2*i,2*j-1)=v*b(m,i)*c(m,j)+(1-v)*c(m,i)*b(m,j)/2;Ke(2*i,2*j)=c(m,i)*c(m,j)+(1-v)*b(m,i)*b(m,j)/2;endendKe=E*t/(4*(1-v^2)*A(EN)).*Ke; %获得单元m的刚度矩阵Kb=mat2cell(Ke,ones(1,3)*2,ones(1,3)*2); %将单元矩阵Ke分为3*3块set1=ones(1,NN)*2;Ka=mat2cell(zeros(NN*2,NN*2),set1,set1); %将总刚K分为NN*NN块for i=1:1:3for j=1:1:3Ka(ecode(m,i),ecode(m,j))=Kb(i,j); %各单元刚度矩阵整体编号,并叠加endendK=K+cell2mat(Ka);end%分块矩阵K合成一个矩阵K%% 引入位移向量和右端项function x=xf(NN,RN,RC,PN,PC,K)x=ones(NN*2,1); %位移初始为0向量for i=1:RNx(RC(i)*2-1)=0;x(RC(i)*2)=0;end %被约束结点位移为0%%----------------引入已知结点载荷-------------%px=zeros(NN*2,1); %载荷初始为0向量for i=1:PNif PC(3,i)==1px(PC(1,i)*2)=PC(2,i);else if PC(3,i)==0px(PC(1,i)*2-1)=PC(2,i);endendend%%----------------引入已知结点载荷-------------%%%----------------求解非约束结点的位移X-------------%set1=ones(1,NN)*2;Ka=mat2cell(K,set1,set1);pxa=mat2cell(px,set1,1);AN=zeros(2*(NN-RN),2*(NN-RN));ANa=mat2cell(AN,ones(1,NN-RN)*2,ones(1,NN-RN)*2);bn=zeros(2*(NN-RN),1);bna=mat2cell(bn,ones(1,NN-RN)*2,1);BN=zeros(2*RN,2*(NN-RN));BNa=mat2cell(BN,ones(1,RN)*2,ones(1,NN-RN)*2);m=1;for i=1:1:NNif x(2*i)==1M(m)=i;m=m+1;endendfor i=1:1:NN-RNfor j=1:1:NN-RNANa(i,j)=Ka(M(i),M(j));bna(i,1)=pxa(M(i),1);endendfor i=1:RNfor j=1:NN-RNBNa(i,j)=Ka(RC(i),M(j));endendAN=cell2mat(ANa);bn=cell2mat(bna);BN=cell2mat(BNa);X=Jordan(AN,bn); %利用Gauss-Jordan法求解非约束结点的位移X %%----------------求解非约束结点的位移X-------------%%----------------由X可得所有结点位移x-------------%BN=BN*X;m=1; n=1;for i=1:1:NNif x(2*i)==1x(2*i-1)=X(m);x(2*i)=X(m+1);m=m+2;else if x(2*i)==0px(2*i-1)=BN(n);px(2*i)=BN(n+1);n=n+2;endendend%% 列主元Jordan消去法将系数矩阵化成对角矩阵求解方程组的数值解function x=Jordan(A,b)%开始计算，赋初值[n,m]=size(A);x=zeros(n,1);for k=1:n%选主元max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endend%交换两行if r>kfor j=k:nz=A(k,j); A(k,j)=A(r,j); A(r,j)=z;endz=b(k); b(k)=b(r); b(r)=z;end%消元计算b(k)=b(k)/A(k,k);for j=k+1:nA(k,j)=A(k,j)/A(k,k);endfor i=1:nif i~=kfor j=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend%输出xfor i=1:nx(i)=b(i);end%% 求解应力应变function [strain stress]=ss(E,v,EN,ecode,A,b,c,x)strain=zeros(3,EN); %应变初始为0矩阵stress=zeros(3,EN); %应力初始为0矩阵D=E/(1-v^2)*[1 v 0;v 1 0;0 0 (1-v)/2];for m=1:1:ENB=[b(m,1) 0 b(m,2) 0 b(m,3) 0;0 c(m,1) 0 c(m,2) 0 c(m,3);c(m,1) b(m,1) c(m,2) b(m,2) c(m,3) b(m,3)];B=B/2/A(m,1); %应变矩阵S=D*B; %应力矩阵X=[x(2*ecode(m,1)-1),x(2*ecode(m,1)),x(2*ecode(m,2)-1),x(2*ecode(m,2)),x(2*ecode(m,3)-1),x (2*ecode(m,3))]';strain(:,m)=B*X;stress(:,m)=S*X;end。

简支梁的有限单元法分析-三角形三节点单元

简支梁的有限单元法分析
三角形三节点平面单元
王峰
有限元分析的基本步骤：
结构离散化
单元分析
整体分析
1 结构离散化
图示为简支梁，梁的厚度为t，泊松比m =0.3，弹性模量为E=2e+5Mpa，用三节点三角形单元进行离散，直角三角形边长为2。
2 单元分析
单元分析的主要内容：由节点位移求内部任一点的
物理方程
{s }=[D]{} 而 { }=[B]{}e (求应力的表达式) {s }=[D][B]{ }e
记 [S]=[D][B]
[S]应力矩阵: [S]=[Si Sj Sm]
2.5节点力与节点位移的关系
令实际受力状态在虚位移状态上做虚功，虚功方程：
({ *}e )T {F}e { *}T tdxdy s
位移，由节点位移求单元应变、应力和节点力。
单元分析的步骤：
节点 (1) 位移
单元内部各点位移
(2)
单元 (3) 应变
单元应力
(4)
节点力
单元分析
2.1 形函数
形函数反映了单元的位移形态，是坐标的函数。三节点三角形单元的形函数为：
1 Ni ( x, y ) (ai bi x ci y )（i , j , m） 2A ai x j ym xm y j （i , j , m） bi y j ym ci xm x j
Ni 1 Ni 1 bi , ci x 2 A y 2 A
因此，三角形单元的应变矩阵[B]是常量，
（i , j , m）
代入数据得到：
1 0 0 0 1 0 1 B 0 0 0 1 0 1 2 0 1 1 0 1 1

有限元三角形单元

有限元三角形单元
有限元三角形单元是有限元分析中常用的基本构件，用于建立复杂结构的数学模型，以进行工程分析。

这些三角形单元是用来离散化连续问题的，将其转换为有限元模型，以便计算机可以进行数值求解。

在有限元分析中，三角形单元有不同类型，其中常见的包括：
1. 线性三角形单元：
- 由三个节点组成的简单三角形单元。

- 三角形的三条边都是直线。

- 这种单元的形状函数是线性的，适用于简单的结构和问题。

2. 二阶和高阶三角形单元：
- 包括更多节点以提高精度的三角形单元。

- 二阶三角形单元具有额外的中间节点，使得形状函数更复杂，从而提高了精度。

- 高阶三角形单元同样通过增加节点来提高精度，但也增加了计算复杂度。

这些三角形单元被用于建立有限元模型，将结构或物体分割成小的几何形状，每个形状都有对应的节点和单元连接关系。

通过这些节点之间的位移和边界条件来建立结构的数学模型，从而进行力学、热力学等各种工程分析。

选择适当类型和精度的三角形单元对于准确地模拟实际问题非常重要，因为它们直接影响到模型的计算精度和效率。

三节点三角形单元算例

三节点三角形单元算例1.引言1.1 概述概述部分（1.1）概述部分旨在介绍本篇文章的主题和背景信息。

本文将重点讨论三节点三角形单元算例，并对其进行详细的分析和结果总结。

三节点三角形单元是有限元分析中常用的一个基本单元，用于建模和模拟复杂的结构系统。

它由三个节点和三条连接这些节点的边组成。

通过这些节点和边的组合，我们可以将结构系统离散成数学模型，进而进行计算和分析。

本文将首先介绍节点的定义和作用。

节点是结构系统中的一个关键概念，它代表了结构的局部特点和重要信息。

在有限元分析中，节点不仅仅是一个几何点，更是一个存储了与该节点相关信息的数据点。

节点的定义和位置决定了单元的布局以及整个结构系统的建模精度。

接着，本文将详细讨论三角形单元的定义和特征。

三角形单元是由三个节点和三条连接这些节点的边组成的简单形状。

这种单元具有较好的适应性和计算效率，广泛应用于各类结构系统的模拟和分析。

通过对三角形单元的了解，我们可以更好地理解和应用有限元方法进行结构分析。

在本文的结论部分，将以一个具体的算例为例，展示三节点三角形单元的应用和分析过程。

通过该算例的分析，我们将得出一些结论和总结，并对三节点三角形单元的适用性和精度进行评估。

在本文的后续章节中，将对节点定义、三角形单元定义、算例分析以及结果总结等内容进行详细阐述，并对相关问题进行讨论和分析。

通过全面的介绍和讨论，我们旨在提供一个全面、准确和有用的指南，帮助读者更好地理解和应用三节点三角形单元算例。

文章结构部分的内容可以按如下方式编写：1.2 文章结构本文共分为引言、正文和结论三个部分。

引言部分首先概述了本文的背景和意义。

随后，介绍了文章的结构，包括本文主要从节点定义和三角形单元定义两个方面展开，以及结语部分。

正文部分主要包括节点定义和三角形单元定义两个子部分。

节点定义部分介绍了在三节点三角形单元中节点的定义和作用。

三角形单元定义部分则详细介绍了三角形单元的结构和数学表示，以及其在有限元分析中的应用。

matlab 有限元基础

Matlab 有限元基础什么是有限元法有限元法（Finite Element Method，FEM）是一种数值计算方法，常被用于工程和科学领域中的结构力学、流体力学、热传导等问题的求解。

有限元法通过将复杂的实际问题离散化为有限个简单的单元，利用数学模型和计算方法来近似求解问题。

有限元法的基本思想是将计算域划分为有限个小单元，每个小单元的物理性质通过节点上的数学函数进行近似描述。

通过对这些小单元的数学模型进行积分计算，得到整个计算域的方程，并通过求解这些方程来得到问题的近似解。

Matlab 在有限元分析中的应用Matlab作为一种功能强大的数值计算和编程软件，被广泛应用于有限元分析中。

它提供了丰富的数学和计算工具，能够方便地实现有限元法的建模、求解和分析。

有限元法的建模在Matlab中，有限元法的建模主要包括以下几个步骤：1.创建几何模型：通过定义节点和单元来描述计算域的几何形状。

可以使用Matlab提供的图形界面工具或者编程方式来创建几何模型。

2.定义边界条件：根据实际问题的边界条件，为模型的节点或单元指定相应的约束条件。

这些条件通常包括位移、力和温度等。

3.定义材料性质：根据实际问题的材料性质，为模型的节点或单元定义相应的材料参数。

这些参数包括弹性模量、泊松比和热导率等。

4.网格划分：将计算域划分为有限个小单元，形成离散化的网格结构。

在Matlab中，可以使用自带的网格划分工具或者自定义的算法进行网格划分。

5.构建刚度矩阵和载荷向量：根据节点的约束条件和单元的材料性质，利用数学公式和计算方法构建刚度矩阵和载荷向量。

有限元法的求解有限元法的求解主要包括以下几个步骤：1.组装系统方程：根据刚度矩阵和载荷向量，将节点的位移和载荷进行组合，形成整个系统的方程。

2.施加边界条件：将已知的位移和载荷应用于系统方程中的相应位置，形成含有未知位移的方程。

3.求解方程：使用Matlab提供的线性代数求解函数，求解含有未知位移的系统方程，得到位移的近似解。

有限元的matlab编程

考虑几何非线性。本程序采用荷载分级加载的方式考虑网架的几何非线性。将总荷载分成1000份分步施加，求解各荷载步下的节点位移，修改网架相应节点坐标以及刚度矩阵，依次迭代求出网架的总位移。
本程序的网架位移求解函数附在主程序后面，主程序运行时调用该函数。
完整ppt
14
几何建模定义荷载
用自定义输入
加下划线的为单元编号
完整ppt
5
形成等效荷载列阵
f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度，a 为之前输入的节点力
集成总刚：
获取单元两端节点坐标
xi = Node( Element( ie, 1 ), 1 ) ;%ie为单元号，以下相同 yi = Node( Element( ie, 1 ), 2 ) ; xj = Node( Element( ie, 2 ), 1 ) ; yj = Node( Element( ie, 2 ), 2 ) ;
节点号：',num2str(Element(ie,2)),' 轴力:',num2str(nodef(1))] ) ;
完整ppt
12
end
例二：网架
完整ppt
13
思路分析
几点说明
网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。鉴于网架的形式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供一种正放四角锥网架的简易输入方法作为典型。
有限元编程示例
完整ppt
1
例一：桁架
题目描述:
如下图所示的平面桁架，杆件长度、弹性模量、截面积以及所受节点力P的大小可以自行定义。求节点位移及杆件轴力。

三角形3节点单元刚度矩阵的具体计算

三角形3节点单元刚度矩阵的具体计算首先，我们需要了解三角形3节点单元的几何形状及节点编号。

三角形3节点单元有三个节点，每个节点有两个坐标值，分别是(x1,y1)，(x2,y2)，(x3,y3)。

对应的节点编号依次为1，2，3接下来，我们将通过数学推导的方式得到三角形3节点单元的刚度矩阵。

1.定义本地坐标系为了简化计算，我们首先引入本地坐标系。

本地坐标系的原点位于三角形的重心，x轴正方向与第1节点连线的方向一致，y轴正方向与x轴正方向在平面内法向量方向一致。

在本地坐标系下，我们可以将三角形3节点单元表示为：```x1=0,y1=0x2=L,y2=0x3=a,y3=b```其中，L是第1节点到第2节点的距离，a和b是第3节点相对于第1节点的坐标。

2.计算本地坐标系下的刚度矩阵在本地坐标系下，我们可以得到单元刚度矩阵Ke的表达式：```Ke=λ*[(b^2+L^2)/4,-a*b/4,-(b^2+L^2)/4,a*b/4,(b^2+L^2)/4,-a*b/4;-a*b/4,(a^2+L^2)/4,a*b/4,-(a^2+L^2)/4,-a*b/4,(a^2+L^2)/4;-(b^2+L^2)/4,a*b/4,(b^2+L^2)/4,-a*b/4,-(b^2+L^2)/4,a*b/4;a*b/4,-(a^2+L^2)/4,-a*b/4,(a^2+L^2)/4,a*b/4,-(a^2+L^2)/4;(b^2+L^2)/4,-a*b/4,-(b^2+L^2)/4,a*b/4,(b^2+L^2)/4,-a*b/4;-a*b/4,(a^2+L^2)/4,a*b/4,-(a^2+L^2)/4,-a*b/4,(a^2+L^2)/4]```其中，λ是三角形的面积，可以通过海伦公式计算：```s=(L+a+b)/2λ = sqrt(s * (s - L) * (s - a) * (s - b))```3.从本地坐标系转换到全局坐标系最后，我们需要将本地坐标系下的刚度矩阵转换到全局坐标系下。

【MATLAB算例】基于3节点三角形单元的矩形薄板分析

【MATLAB 算例】基于 3节点三角形单元的矩形薄板分析将此结构按三角形单元划分成 432个三角形（X 方向分成 18段， Y 方向分成12 段），总共分成 19X13=247个结点的有限元模型，具体步骤详细程序如下：tic;Initial_info=[0.09 0.06 18 12];disp(该程序计算的是 ',num2str(Initial_info(3)+1),'X',num2str(Initial_info(4)+1),'=',...num2str((Initial_info(3)+1)*(Initial_info(4)+1)),' 个结点的有限元模型 ']);LX=Initial_info(1);LY=Initial_info(2);nx=Initial_info(3);ny=Initial_info(4);ne=2*nx*ny;np=(nx+1)*(ny+1);for i=1:nx+1;j=1:ny+1;Np(i,j)=j+(i-1)*(ny+1);end生成节点编号矩阵 Npfor i=1:nx+1;j=1:ny+1;XX(i,j)=(i-1)*LX/nx;YY(i,j)=(j-1)*LY/ny;endXY=[reshape(XX',np,1),reshape(YY',np,1)];nx2=nx/2;Np1=Np(1:nx2+1,:);Np2=Np(nx2+1:end,:);for i=1:nx2*ny;if rem(i,nx2)==0xp=nx2;yp=i/nx2;elsexp=rem(i,nx2);yp=fix(i/nx2)+1;endDof1(i,:)=[Np1(xp,yp),Np1(xp+1,yp),Np1(xp,yp+1)];Dof1(i+nx2*ny,:)=[Np1(xp+1,yp),Np1(xp+1,yp+1),Np1(xp,yp+1)];Dof2(i,:)=[Np2(xp,yp),Np2(xp+1,yp),Np2(xp+1,yp+1)];Dof2(i+nx2*ny,:)=[Np2(xp,yp),Np2(xp+1,yp+1),Np2(xp,yp+1)];endDof=[Dof1;Dof2];for i=1:neunit(i,:)=[XY(Dof(i,1),1),XY(Dof(i,2),1),XY(Dof(i,3),1),...XY(Dof(i,1),2),XY(Dof(i,2),2),XY(Dof(i,3),2)];enddisp('前处理完成 ');前处理完成单元刚度矩阵E=2*10^11;u=0.3;平面应力问题D=E/(1-u^2)*[1 u 0;u 1 0;0 0 (1-u)/2];for i=1:nexi=unit(i,1);yi=unit(i,4);xj=unit(i,2);yj=unit(i,5);xm=unit(i,3);ym=unit(i,6);ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi;bi=yj-ym;bj=ym-yi;bm=yi-yj;ci=-(xj-xm);cj=-(xm-xi);cm=-(xi-xj);area=abs((ai+aj+am)/2);B = [bi 0 bj 0 bm00 ci 0 cj0 cmci bi cj bj cm bm];Be{i,1} = B/2/area ;ke{i,1}=[Be{i,1}]'*D*Be{i,1}*area;end总刚度矩阵叠加KK=sparse(2*np,2*np);for ie=1:nea=Dof(ie,1);b=Dof(ie,2);c=Dof(ie,3);DOF(1)=2*a-1;DOF(2)=2*a;DOF(3)=2*b-1;DOF(4)=2*b;DOF(5)=2*c-1;DOF(6)=2*c;for n1=1:6for n2=1:6KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+ke{ie,1}(n1,n2);endendend单元等效节点荷载y=(0:LY);P=(10^7/0.03)*y-10^7; 左右受变化的三角形荷载 ,在如图的坐标系下Re=sparse(ne,6);for i=1:ne;switch icase num2cell(1:ne/2-nx2)Pe=[0 0 0 0 0 0];case num2cell(ne/2-nx2+1:ne/2)Pe=-LX*P*[0,0,0,1,0,1]/nx/2;case num2cell(ne/2+1:ne-nx2)Pe=[0 0 0 0 0 0];otherwisePe=-LX*P*[0,0,0,1,0,1]/nx/2;endRe(i,:)=Pe;end荷载叠加Rr=sparse(1,2*np);for i=1:nea=Dof(i,1);b=Dof(i,2);c=Dof(i,3);DOF(1)=2*a-1;DOF(2)=2*a;DOF(3)=2*b-1;DOF(4)=2*b;DOF(5)=2*c-1;DOF(6)=2*c;for n1=1:6Rr(DOF(n1))= Rr(DOF(n1))+Re(i,n1);endend生成需处理的行列cp=[1:nx+1];ctype=ones(1,length(cp));ctype(nx2+1)=2;cp_all=(cp-1)*(ny+1)+1;p_stake=zeros(1,2*length(cp));for i=1:length(cp)switch ctype(i)case {2}p_stake(2*i)=2*cp_all(i);p_stake(2*i-1)=2*cp_all(i)-1;case {1}p_stake(2*i)=2*cp_all(i);p_stake(2*i-1)=[];otherwisep_stake(2*i)=[];p_stake(2*i-1)=2*cp_all(i)-1;endend[m,j]=find(p_stake==0);p_stake(:,j)=[];处理对应的行列KK_d=KK;KK_f=KK;KK_d(p_stake,:)=[];KK_d(:,p_stake)=[];KK_f(:,p_stake)=[];KK_f=KK_f(p_stake,:);Rr_unkown=Rr;Rr_unkown(:,p_stake)=[];RR=transpose(Rr_unkown);[L,U]=lu(KK_d);UU=U\(L\RR);Rx=KK_f*UU;数值计算部分UU_all=UU';for i=1:length(p_stake)UU_all=[UU_all(:,1:p_stake(i)-1),0,UU_all(:,p_stake(i):end)];endfor i=1:npUU_info(i,:)=[UU_all(2*i-1),UU_all(2*i)];end出图部分（运行后显示）梁尺寸及荷载图figure;set(gcf,'outerposition',get(0,'ScreenSize'));set(gcf,'name','梁的尺寸及荷载 ');line([-0.015,-0.01],[0.06,0.06]),hold online([-0.01,-0.01],[0.06,0]),hold online([-0.01,-0.005],[0,0]),hold online([-0.005,-0.015],[0,0.06]),hold online([0.105,0.1],[0.06,0.06]),hold online([0.1,0.1],[0.06,0]),hold online([0.1,0.095],[0,0]),hold online([0.095,0.105],[0,0.06]),hold onrectangle('position',[0,0,0.09,0.06]),hold on quiver(-0.01,0.06,-0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.05,-0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold onquiver(-0.01,0.04,-0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.02,0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.01,0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0,0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.06,0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.05,0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.04,0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.02,-0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.01,-0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0,-0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold ontext(-0.0155,0.062,'1000N/cm^2')text(-0.0105,-0.002,'1000N/cm^2')text(0.1,0.062,'1000N/cm^2')text(0.095,-0.002,'1000N/cm^2')text(0,-0.002,'梁宽 9cm，高 6cm，厚 1cm。

《有限元基础教程》_【MATLAB算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)

如三节点三角形单元

《有限元分析》课程作业

有限元程序设计--第五章 线性三角形单元

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

基于MTLAB的3节点三角形单元求解平面弹性问题

matlab有限元三单元编程

有限元分析基础教程

matlab有限元三角形单元编程

matlab有限元三角形单元程序

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

《有限元基础教程》_【MATLAB算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(Triangle2D3Node)

3结点三角形单元有限元程序MATLAB语言

简支梁的有限单元法分析-三角形三节点单元

有限元三角形单元

三节点三角形单元算例

matlab 有限元基础

有限元的matlab编程

三角形3节点单元刚度矩阵的具体计算

【MATLAB算例】基于3节点三角形单元的矩形薄板分析

有限元程序设计--第五章线性三角形单元