平面三角形单元有限元程序设计.docx

1、平面三角形单元有限元程序设计、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖 向载荷。已知：P=150N/m , E=200GPa ，=0.25 , t=0.1m ，忽略自重。试计算薄板的位移及应力分布。要求:1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）；2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程 序计算结果进行对比（任选取三个点，对比位移值）；3.提交程序编写过程的详细报告及计算机程序；4.所有同学参加答辩，并演示有限元计算程序。有限元法

2、中三节点三角形分析结构的步骤如下：1） 整理原始数据，如材料性质、荷载条件、约束条件等，离散结构 并进行单元编码、结点编码、结点位移编码、选取坐标系。2） 单元分析，建立单元刚度矩阵。3） 整体分析，建立总刚矩阵。4） 建立整体结构的等效节点荷载和总荷载矩阵5） 边界条件处理。6） 解方程，求出节点位移7） 求出各单元的单元应力8） 计算结果整理。、程序设计网格划分如图，将薄板如图划分为6行，并建立坐标系，则节点编号单兀编号刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个 单元对应的节点序号。通过循环逐个计算：（

3、1 ）每个单元对应2种单元刚度矩阵中的哪一种；（2） 该单元对应总刚度矩阵的那几行哪几列（3） 将该单元的单元刚度矩阵加入总刚度矩 阵的对应行列循环又分为3层循环：（1）最外层：逐行计算（2） 中间层：该行逐个计算（3） 最里层：区分为第奇/偶数个计算k1e66k1 56 56k2e66k2 56 56单元刚度的集成： kZ66kZ 56 56Kke k e12kZe边界约束的处理：划 0 置 1 法适用： 这种 方法 适用 于 边界 节点 位移 分量 为已 知(含 为 0) 的各 种约 束。做法：(1 )将总刚矩阵K中相应于已知位移行主对角线元素置1 ,其 他元素改为零；同 时将载荷列阵 R

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

5、 总 刚 元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单元应力%*%初始化clearformat short e%设定输 出类型clear%清

6、除内存变量NELEM=36%单 元 个 数 ( 单 元 编 码 总 数 )NPION=28%结 点 个 数 ( 结 点 编 码 总 数 )NVFIX=2%受 约 束 边 界 点 数NFORCE=1%结点荷载个数YOUNG=2e11%弹性模量POISS=0.25%泊松比%厚度THICK=0.1LNODS=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; 12 18 13;

7、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.5 6;3 6

8、;-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，无约%*%生 成单 元 刚度矩阵 并组成总 体刚度矩 阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总 体刚 度矩阵 并置 0

9、%*for i=1:NELEM%生 成弹 性 矩阵 DD= 1POISS 0;POISS1 0;00 (1-P0ISS)/2*Y0UNG心-POISST)%*%* %生 成应 变 矩阵 Bfor j=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);endB=b(1)0b(2)0b(3)0;0c(1)0c(2)0c(3);c(1)b(1) c(2

10、)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%*%将约

11、束信息加入总体刚度矩阵(对角元素改一法)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;% 对 角 元素 为 1end%*%生 成单 元 刚度矩阵 并组成总 体刚度矩 阵%* if FIXED(i,3)=1ASTIF( :,FIXED(i,1)*2)=0; %一列 为零ASTIF(FIXED(i,1)*2,:)=0; %一行 为零ASTIF(FIXED(i,1)*2

12、 ,FIXED(i,1)*2)=1; % 对 角 元 素 为1endend%*%生成荷载向量ASLOD(1:2*NPION)=0; %总体 荷载向 量置零for i=1:NFORCE ASLOD(FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%*%求 解内 力ASDISP=ASTIFASLOD %计算 节点 位 移向 量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);%取出当前单元

13、的节点位移向量endi%求内力STRESS=D*B1(:, :, i)*ELEDISPend（程序计算结果和有限元软件得出的结果稍有偏差，可能是程序某些 地方数据输入时出了问题，还在寻找具体原因）有限元软件分析1PropertyhaliteUnit5铠| Referencie Temperature226- 诃 Isotropic Elasbdty7Derive fromYoirigs mi T8Youngs Modulus2E+11P3 工gPatssons Ratio0.2510BiA Modulusr 1.3333E+11Pa11Shear Modulusfif+iOPa12Fl 泊 F

14、ield Variables-设置材料参数U&tflils ot Su rtarff BodyGraphics Prcperti-esA-DefinitionSuppressedNOStiffnen BehaviorFlexibleCoord mate SystemMau It Coordinate SystemReference TemperatureBy EnvironmentThicknessIGO. mmThickness MtdeManualOffset TypwMiddle-MaterialAwignmntStrurturaJ StwlPJortlinear E ffectsYes

15、V SelectiorGeomettyiBodyDefiniitiariSdpprsdNoMethodTriangles|EenientMid de INodesUse Global Sdtn j 寸Z)etai s o: All Triangles Method MethodGecmrtry3Eog-DefhitionSupprwsedhkjrnMrofDivsicnsNdNunbo2 DvsionsE#c7I0T-ardEdNd Biasc_ .斤rm网格划分雪 一 2 曾(wnj边界约束添加载荷A： Katie ttractural匚qjivdlcri. Sb 4札-qL.ii/ pm

16、卜GlvMirpc Strns - Top./Rntrnin Unite MPc-irrifl： I2flT7/ip 11;14CQJ0001 S91 M MaxQX0C1&93； tjNCUMDOJQOCl?DJJOCIMfl?099E 讥 55JO32C-30jC4G7eW6d SMiftA: Malic StnictunlTci柄| 口 frfoFTTiflranType; -Q7fll DF:jrQi*i0iUnb*: mmTime:3017/1/7 11:12Jt6K1e-6 Mai4.1716c &cxb5D2e-63.1 237e-62XD73e52jOB5Sb-&k5W4e-61.04?9e-&5L2145b-70 Min