有限元法MATLAB作业

例题：
如图所示的受力均匀分布载荷作用的薄平板结构，将平板离散成两个线性三角元，如图所示，假定2
210,0.3,0.025,3000/m E GPa v t m w KN ====。

求 （1）该结构的整体刚度矩阵。

（2）节点2和节点3的水平位移和垂直位移。

（3）节点1和节点4的支反力。

（4）每个单元的应力。

（5）每个单元的主应力和主应力方向角。

图2. 用双线性三角元离散化的薄木板 解：（1）离散化域
我们将平板分为两个单元，4个节点，如图2所示，分布载荷的总作用力平均分给节点2和节点3，由于结构是薄平板，所以假定其属于平面应力情况。

MATLAB 中采用的计算单位是KN 和m 。

表1给出了该题的单元连通性。

（2）单元刚度矩阵
通过调用MATLAB 的LinearTriangleElementStiffness 函数，得到两个单元矩阵k 1和k 2，每个矩阵都是66⨯的。

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(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)
图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 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.2740
k2=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,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 -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）集成整体刚度矩阵
⨯的，因此为了得到整体刚度矩阵K，我们由于结构有4个节点，所以刚度矩阵是88
⨯的零矩阵。

由于在这个系统中有两个线性三角元，所以两次调用MATLAB 首先要生成一个88
的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 0 0 0 0 0 0
K1 = 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
0 5.7692 0 0 -0.8654 0 0.8654 -5.7692
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 -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.8750
1.0096 -5.7692 0 0 0.8654 -0.5048 -1.8750 6.2740 K2 = gdfuhe2( @ LinearTriangleAssemble,K1,k2,1,2,3 )的源程序：
function 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.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）引入边界条件
用上一步得到的整体刚度矩阵，可以得到该结构的方程组，如下所示
6 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.00910 6 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 111122233441.0096 3.4615 -1.8750 1.0096 -5.7692 0 0 0.8654 -0.5048 -1.8750 6.2740x x y y x y x y x y U F U F U F U U U U U ⎧⎫⎡⎤⎪⎪⎢⎥⎪
⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥=⎨⎬⎢⎥⎪
⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎣⎦⎩⎭23344x y x y x y F F F F F ⎧⎫⎪⎪
⎪⎪⎪⎪⎪⎪
⎪⎪⎨⎬
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭
（1）
本题的边界条件如下：
114422330
9.375,0,9.375,0
x y x y x y x y U U U U F F F F ========
将上述边界条件代入方程（1）中，可以得到：
6 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.0096 10 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 311223344009.37509.3750
.4615 -1.87500 1.0096 -5.7692 0 0 0.8654 -0.5048 -1.8750 6.27400x y x y x
y x y F F U U U U F
F ⎧⎫⎧⎡⎤⎪⎪⎪⎢⎥⎪
⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥=⎨⎬⎨⎢⎥⎪
⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎪⎪⎢⎥⎣⎦⎩⎭⎫
⎪⎪⎪⎪
⎪
⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭
（2）
（5）解方程及后处理
用分解（手动）和高斯消去法（MATLAB ）求解方程组。

首先对方程组（2）进行分解，提取整体刚度矩阵K 的第3~6行和第3~6列作为子矩阵。

从而得到：
22633 3.4615 -1.8750 -2.0192 0.8654 1.8750 6.2740 1.0096 -5.769210-2.0192 1.0096 3.4615 0 0.8654 -5.7692 0 6.2740x y x y U U U U ⎧⎫⎡⎤⎪⎪⎢⎥⎪⎪⎢⎥=⎨⎬⎢⎥⎪⎪⎢⎥⎪⎪⎣⎦⎩⎭9.37509.3750
⎧⎫⎪⎪
⎪⎪⎨⎬
⎪⎪⎪⎪⎩⎭ 方程组的解用如下MATLAB 命令可以得到。

[u,U,F]=weiyiyingli( @ LinearTriangleAssemble,K1,k2,1,2,3 )的源程序：
function [u,U,F]=weiyiyingli( f,K,k,i,j, m )
%UNTITLED3 此处显示有关此函数的摘要
% 此处显示详细说明[u,U,F]=weiyiyingli( @ LinearTriangleAssemble,K1,k2,1,2,3 )
K2=f(K,k,i,j,m);
k=K2(3:6,3:6);
f=[9.375;0;9.375;0];
u=k\f;
U=[0;0;u;0;0];
F=K2*U;
>> [u,U,F]=weiyiyingli( @ LinearTriangleAssemble,K1,k2,1,2,3 )
u =
1.0e-05 *
0.7111
0.1115
0.6531
0.0045
U =
1.0e-05 *
0.7111
0.1115
0.6531
0.0045
F =
-9.3750
-5.6295
9.3750
0.0000
9.3750
0.0000
-9.3750
5.6295
从u中可得：节点2的水平位移和垂直位移分别是0.7111m和0.1115m。

节点3的水平位移和垂直位移分别是0.6531m和0.0045m。

从F中可得：节点1的水平支反力和垂直支反力分别是9.375KN（指向左边）和5.6295KN （作用力方向向下），节点4的水平支反力和垂直支反力分别是9.375KN（指向左边）和5.6295KN（作用力方向向上），满足力的平衡。

接着，我们建立单元节点位移矢量u1和u2，然后调用MATLAB的LinearTriangleElementStresses函数计算单元的应力sigmal1和sigmal2,及调用MA TLABA的LinearTriangleElementPStresses函数计算每个单元的主应力和主应力方向角。

u1u2的源程序：
u1=[U(1);U(2);U(5);U(6);U(7);U(8)]
u2=[U(1);U(2);U(3);U(4);U(5);U(6)]
sigma1=LinearTriangleElementStresses(210000000, 0.3, 0, 0, 0.5, 0.25, 0, 0.25, 1, u1)
sigma2=LinearTriangleElementStresses(210000000, 0.3, 0, 0, 0.5, 0 , 0.5, 0.25, 1, u2)
s1 = LinearTriangleElementPStresses(sigma1)
s2 = LinearTriangleElementPStresses(sigma2)
>> u1u2
u1 =
1.0e-05 *
0.6531
0.0045
u2 =
1.0e-05 *
0.7111
0.1115
0.6531
0.0045
sigma1 =
1.0e+03 *
3.0144
0.9043 0.0072
sigma2 =
1.0e+03 *
2.9856 -0.0036 -0.0072 s1 =
1.0e+03 *
3.0144 0.9043 0.0002 s2 =
1.0e+03 *
2.9856 -0.0036 -0.0001
从sigma1和sigma2中可得：单元1的应力 3.0144x MPa σ=（拉应力），
0.9043y MPa σ=（拉应力），0.0072xy MPa τ= （正值）。

单元1的应力 2.9856x MPa σ=（拉应力），
0.0036y MPa σ=（压应力），0.0072xy MPa τ=（负值）。

很显然，在x 方向的应力（拉
应力）接近于正确的值3MPa （拉应力）。

从s1和s2中可得：单元1的应力 3.0144x MPa σ=（拉应力），0.9043y MPa σ=（拉应力），主应力方向角0.2p θ=︒ 。

单元1的应力 2.9856x MPa σ=（拉应力），
0.0036y MPa σ=（压应力），主应力方向角0.1p θ=-︒。