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

基于MTLAB的3节点三角形单元求解平面弹性问题 引言 本程序相比其他参考资料的程序，在操作方面得到大大简化，只需要一次输入相关力学参量、单元总数、节点总数、节点坐标、单元节点编号，利用此程序就能够很快得到所有单元的单元刚度矩阵，以及系统的整体刚度矩阵。为了便于后面计算单元的应力，在计算单元刚度矩阵时，MATLAB程序输出了弹性系数矩阵D和几何矩阵B。 此外，本程序在引入位移边界条件求解刚度方程时，采用了化零置1法，避免了手工化简整体刚度矩阵的麻烦。当然，这在一定程度上加重了计算机的负担，增加了计算机的计算时间，但是这个延长的时间相对于手工化简刚度矩阵所需的时间来说，是微不足道的。 程序需要提前输入的参数有： 节点坐标node_gcoord 单元节点编号element_node_number 节点总数node_number 单元总数element_number 弹性模量E 泊松比NU 厚度t 平面问题性质指示参数ID，其中1为平面应力问题，2为平面应变问题 该程序主要包括以下部分： 单元刚度矩阵程序：用于计算各单元的单元刚度矩阵； 整体刚度矩阵程序：用于组装整体刚度矩阵； 位移及支反力求解程序：用于求解未知位移和未知支反力 单元应力计算程序：用于计算各单元的单元应力 问题描述 为了验证程序的正确性，我们选取了参考文献[1]中的一个简单实例进行验证。问题如下： 如图所示为一矩形薄平板，在右端部受集中力F=100kN的作用，材料常数为：弹性模量E=1e7Pa，泊松比NU=1/3，板的厚度t=0.1m。试按平面应力问题计算各节点位移及支座反力。 问题求解 结构的离散化与编号 节点描述 2) 单元描述 计算单元刚度矩阵 KE(:,:,1) = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.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 KE(:,:,2) = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.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 组装整体刚度矩阵 KK = 1.0e+006 * 0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 0 0 0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 0 0 -0.3750 -0.1875 0.6563 0 0 0.3750 -0.2813 -0.1875 -0.1875 -1.1250 0 1.2188 0.3750 0 -0.1875 -0.0938 -0.2813 -0.1875 0 0.3750 0.6563 0 -0.3750 -0.1875 -0.1875