《有限元法MATLAB作业.docx》由会员分享，可在线阅读，更多相关《有限元法MATLAB作业.docx（7页珍藏版）》请在金锄头文库上搜索。
1、例题：如图所示的受力均匀分布载荷作用的薄平板结构，将平板离散成两个线性三角元，如图所示，假定。求（1）该结构的整体刚度矩阵。（2）节点2和节点3的水平位移和垂直位移。（3）节点1和节点4的支反力。（4）每个单元的应力。（5）每个单元的主应力和主应力方向角。图1. 薄平板结构图2. 用双线性三角元离散化的薄木板解：（1）离散化域我们将平板分为两个单元，4个节点，如图2所示，分布载荷的总作用力平均分给节点2和节点3，由于结构是薄平板，所以假定其属于平面应力情况。MATLAB中采用的计算单位是KN和m。表1给出了该题的单元连通性。表1 该题的单元连通性单 元 编 号节 点 i节 点 j节 点 m11
2、342123（2）单元刚度矩阵通过调用MATLAB的LinearTriangleElementStiffness函数，得到两个单元矩阵k1和k2，每个矩阵都是的。k1=gangdujuzhen(LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)的源程序：function k1= gangdujuzhen( f,E,NU,t,xi,yi,xj,yj,xm,ym,p )%UNTITLED4 此处显示有关此函数的摘要% 此处显示详细说明k1=gangdujuzhen(LinearTri
3、angleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)k1= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);k1=gangdujuzhen(LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)k1 = 1.0e+06 * 2.0192 0 0 -1.0096 -2.0192 1.0096 0 5.7692 -0.8654 0 0.8654 -5.7692 0 -0.8654
4、1.4423 0 -1.4423 0.8654 -1.0096 0 0 0.5048 1.0096 -0.5048 -2.0192 0.8654 -1.4423 1.0096 3.4615 -1.8750 1.0096 -5.7692 0.8654 -0.5048 -1.8750 6.2740k2=gangdujuzhen2(LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)的源程序：function k2= gangdujuzhen2( f,E,NU,t,xi,yi,xj,yj,
5、xm,ym,p )%UNTITLED3 此处显示有关此函数的摘要% 此处显示详细说明 k2=gangdujuzhen2(LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);k2=gangdujuzhen2(LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2 = 1.0e+06 * 1.4423 0 -
6、1.4423 0.8654 0 -0.8654 0 0.5048 1.0096 -0.5048 -1.0096 0 -1.4423 1.0096 3.4615 -1.8750 -2.0192 0.8654 0.8654 -0.5048 -1.8750 6.2740 1.0096 -5.7692 0 -1.0096 -2.0192 1.0096 2.0192 0 -0.8654 0 0.8654 -5.7692 0 5.7692（3）集成整体刚度矩阵由于结构有4个节点，所以刚度矩阵是的，因此为了得到整体刚度矩阵K，我们首先要生成一个的零矩阵。由于在这个系统中有两个线性三角元，所以两次调用MATL
7、AB的LinearTriangleAssemble函数就可以得到整体刚度矩阵K。每次对该函数的调用都会集成一个单元。K = zero( 8,8 )的源程序：function K = zero( x,y )%UNTITLED4 此处显示有关此函数的摘要% 此处显示详细说明 K = zero( 8,8 )K=zeros(x,y); K = zero( 8,8 )K = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
8、0 0 0 0 0 0K1 = gdfuhe1( LinearTriangleAssemble,K,k1,1,3,4 )的源程序：function K1 = gdfuhe1( f,K,k,i,j,m )%UNTITLED 此处显示有关此函数的摘要% 此处显示详细说明 K1 = gdfuhe1( LinearTriangleAssemble,K,k1,1,3,4 )K1=f(K,k,i,j,m); K1 = gdfuhe1( LinearTriangleAssemble,K,k1,1,3,4 )K1 = 1.0e+06 *2.0192 0 0 0 0 -1.0096 -2.0192 1.0096
9、0 5.7692 0 0 -0.8654 0 0.8654 -5.76920 0 0 0 0 0 0 00 0 0 0 0 0 0 00 -0.8654 0 0 1.4423 0 -1.4423 0.8654-1.0096 0 0 0 0 0.5048 1.0096 -0.5048-2.0192 0.8654 0 0 -1.4423 1.0096 3.4615 -1.87501.0096 -5.7692 0 0 0.8654 -0.5048 -1.8750 6.2740K2 = gdfuhe2( LinearTriangleAssemble,K1,k2,1,2,3 )的源程序：function
10、 K2= gdfuhe2( f,K,k,i,j,m )%UNTITLED2 此处显示有关此函数的摘要% 此处显示详细说明 K2 = gdfuhe2( LinearTriangleAssemble,K1,k2,1,2,3 )K2=f(K,k,i,j,m); K2 = gdfuhe2( LinearTriangleAssemble,K1,k2,1,2,3 )K2 = 1.0e+06 * 3.4615 0 -1.4423 0.8654 0 -1.8750 -2.0192 1.0096 0 6.2740 1.0096 -0.5048 -1.8750 0 0.8654 -5.7692 -1.4423 1
11、.0096 3.4615 -1.8750 -2.0192 0.8654 0 0 0.8654 -0.5048 -1.8750 6.2740 1.0096 -5.7692 0 0 0 -1.8750 -2.0192 1.0096 3.4615 0 -1.4423 0.8654 -1.8750 0 0.8654 -5.7692 0 6.2740 1.0096 -0.5048 -2.0192 0.8654 0 0 -1.4423 1.0096 3.4615 -1.8750 1.0096 -5.7692 0 0 0.8654 -0.5048 -1.8750 6.2740（4）引入边界条件用上一步得到的整体刚度矩阵，可以得到该结构的方程组，如下所示本题的边界条件如下： 将上述边界条件代入方程（1）中，可以得到：（5）解方程及后处理用分解（手动）和高斯消去法（MATLAB）求解方程组。首先对方程组（2）进行分解，提取整体刚度矩阵K的第36行和第36列作为子矩阵。从而得到：方程组的解用如下MATLAB命令可以得到。u,U,F=weiyiyingli( LinearTriangleAssemble,K1,k2,1,2,3 )的源程序：function u,U,F=weiyiyingli( f,K,k,i,j, m )%UNTITLED3 此处显示有关此函数的摘
