Matlab有限元分析操作基础

Matlab 有限元分析20140226
为了用Matlab 进行有限元分析，首先要学会Matlab 基本操作，还要学会使用Matlab 进行有限元分析的基本操作。

1. 复习：上节课分析了弹簧系统
x 推导了系统刚度矩阵
11221
21200k k k k k k k k -⎡⎤⎢⎥-⎢⎥⎢⎥--+⎣⎦
2. Matlab有限元分析的基本操作
（1）单元划分（选择何种单元，分成多少个单元，标号）（2）构造单元刚度矩阵（列出…）
（3）组装系统刚度矩阵（集成整体刚度矩阵）
（4）引入边界条件（消除冗余方程）
（5）解方程
（6）后处理（扩展计算）
3. Matlab有限元分析实战【实例1】
分析：
步骤一：单元划分
步骤二：构造单元刚度矩阵
>>k1=SpringElementStiffness(100) >>…?
步骤三：构造系统刚度矩阵
a) 分析SpringAssemble库函数function y = SpringAssemble(K,k,i,j)
% This function assembles the element stiffness
% matrix k of the spring with nodes i and j into the % global stiffness matrix K.
% function returns the global stiffness matrix K
% after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
y = K;
b) K是多大矩阵？
今天的系统刚度矩阵是什么？
因为
11
22
1212
k k
k k
k k k k
-
⎡⎤
⎢⎥
-
⎢⎥⎢⎥--+
⎣⎦
所以
1000100
0200200 100200300
-
⎡⎤
⎢⎥
-
⎢⎥
⎢⎥
--
⎣⎦
？
c) K=SpringAssemble(K,k1,1,2) function y = SpringAssemble(K,k,i,j)
K(i,i) = K(i,i) + k(1,1);
K(i,j) = K(i,j) + k(1,2);
K(j,i) = K(j,i) + k(2,1);
K(j,j) = K(j,j) + k(2,2);
1100100100100k -⎡⎤
=⎢⎥-⎣⎦
10010001001000000K -⎡⎤
⎢⎥
=-⎢⎥
⎢⎥⎣⎦
K=SpringAssemble(K,k2,2,3) 1200200200200k -⎡⎤
=⎢⎥-⎣⎦
10010001003002000200200K -⎡⎤
⎢⎥
=--⎢⎥⎢⎥-⎣⎦
1001000100010010030020002002000200200100200300--⎡⎤⎡⎤⎢⎥⎢⎥--≠-⎢⎥⎢
⎥⎢⎥⎢⎥---⎣⎦⎣⎦!?
步骤四：引入边界条件，消除冗余方程>>k=K(2:3,2:3)%构造不含冗余的方程>>f=[0;15]%构造外力列阵
步骤五：解方程
引例：已知1212
u 31u u u +=⎧⎨-=⎩，求 12u u 和 解：
类似求解KU=F ，
输入下列Matlab 命令：
>> K=[1 1;1,-1]
>> F=[3;1]
>> U=inv(K)*F
>> U=K \F
（继续弹簧系统求解）
>>u=k \f %使用高斯消去法求解 >>U=[0 ; u]%构造原方程组
>>F=K*U %求出所有外力，含多余计算
步骤六：后处理、扩展计算
>>u1=[0;U(2)]%构造单元位移
>>f1=SpringElementForces(k1,u1)%求单元1内力>>u2=[U(2) ; U(3)]%构造单元2位移
>>f2=SpringElementForces(k2,u2)%求单元2内力
4. 总结
clc
clear
k1=SpringElementStiffness(100)%创建单元刚度矩阵1 k2=SpringElementStiffness(200)%创建单元刚度矩阵2 K=zeros(3,3)%创建空白整体刚度矩阵
K=SpringAssemble(K,k1,1,2)%按节点装入单元矩阵1 K=SpringAssemble(K,k2,2,3)%按节点装入单元矩阵2 k=K(2:3,2:3)%构造不含冗余的方程
f=[0;15]%构造外力列阵
u=k\f%使用高斯消去法求解
U=[0 ; u]%构造系统节点位移列阵
F=K*U%求出所有外力，含多余计算
u1=[0;U(2)]%构造单元位移
f1=SpringElementForces(k1,u1)%求单元1内力
u2=[U(2) ; U(3)]%构造单元2位移
f2=SpringElementForces(k2,u2)%求单元2内力
5. 练习
1 Danyi 13
2 dan 34 3dan 35 4dan 35 dan5 54 dan6 42。