平面三角形单元有限元程序设计文档格式.doc

1、如图，将薄板如图划分为6行，并建立坐标系，则PXY单元编号节点编号刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算单元刚度的集成：边界约束的处理：划0置1法适用：这种方法适用于边界节点位移分量为已知（含为0）的各种约束。做法：

2、（1）将总刚矩阵K中相应于已知位移行主对角线元素置1，其他元素改为零；同时将载荷列阵R中相应元素用已知位移置换。这样，由该方程求得的此位移值一定等于已知量。（2）将K中已知位移相应的列的非主对角成元素也置0，以保持K的对称性。当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去已知位移对该方程的贡献已知位移和相应总刚元素的乘积。若约束为零位移约束时，此步则可省去。特点：经以上处理同样可以消除刚性位移（约束足够的前提下），去掉未知约束反力。但这种方法不改变方程阶数，利于存贮。（3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。程序如下：变

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 单元应力 %*%初始化clearformat short

4、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;

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

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

7、 for i=1:NELEM %生成弹性矩阵 D D= 1 POISS 0; POISS 1 0; 0 0 （1-POISS）/2*YOUNG/（1-POISS2） %计算当前单元的面积 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（LNO

8、DS（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）; 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

9、=1: ASTIF（a（j）*2-1）:a（j）*2,（a（k）*2-1）:a（k）*2）=ASTIF（a（j）*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