(今天挺忙的先水水吧,代码直接放在这里了,另外附上少量的图解,以后如果有时间我会专门法一篇文章讲有限元)
下面这个是matlab平面四节点等参元的刚度矩阵计算程序
第一个是test.mlx,这个部分是计算刚度矩阵的符号矩阵的,首先对于四节点单元,其编号是逆时针顺序编号的,即:
对于四节点单元,其单元刚度矩阵是8x8的矩阵,其符号矩阵计算如下
test.mlx (这个用于测试, mlx文件只要在左边文件栏点击新建live script就行)
% 四节点单元刚度矩阵函数 E杨氏模量, nu 泊松比, h板厚度
syms x
syms y
syms E nu h
M = [
[1, x, y, x*y,0,0,0,0]
[0,0,0,0, 1, x, y, x*y]
];
A = [
[1 -1 -1 1 0 0 0 0]
[0 0 0 0 1 -1 -1 1]
[1 1 -1 -1 0 0 0 0]
[0 0 0 0 1 1 -1 -1]
[1 1 1 1 0 0 0 0]
[0 0 0 0 1 1 1 1]
[1 -1 1 -1 0 0 0 0]
[0 0 0 0 1 -1 1 -1]
];
N = M/A; % M * inv(A) , 为形状函数矩阵
% 对x, y等参变量求导并得到对应的应变矩阵
% 对x, y等参变量求导并得到对应的应变矩阵
N_x = diff(N,x);
N_y = diff(N,y);
B = [N_x(1,:)
N_y(2,:)
N_y(1,:) + N_x(2,:)];
clear N_x N_y
D = E /(1-nu^2) * [
[1 nu 0]
[nu 1 0]
[0 0 (1- nu)/2]
];
k_e = transpose(B) * D * B;
Ke = h*int(int(k_e,x,[-1, 1]), y, [-1,1]); % 积分求解对应的四节点刚度矩阵
Ke2 = Ke/(E * h/(1-nu^2));
simplify(Ke2)
求解结果是这样的(也可以使用一个.m文件来运行)
这个的正确性可以参考 《结构有限元分析》 刘永寿 第二版 第51页, 上面的刚度矩阵将a/b = 1代入之后和这个刚度矩阵结果相同。
因此我们求解单元刚度矩阵只需要将值代入就行了
下面附上可以直接复制粘贴使用的MATLAB函数代码
这个函数的调用方法是K_e = elem401(E, nu, h)
其中E是杨氏模量,nu是泊松比,h是平面板的厚度
比如elem401(2e5, 0.3, 10)
其参数是杨氏模量200Gpa(由于单位是mm,而Pa是
N
/
m
2
N/m^2
N/m2 ,所以
2
×
1
0
11
P
a
=
2
×
1
0
5
N
/
m
m
2
)
2\times10^{11}Pa=2\times 10^5N/mm^2)
2×1011Pa=2×105N/mm2)
上面的调用方法得到刚度单位的是
N
/
m
m
N/mm
N/mm (你使用mm为单位计算)
当然使用m为单位计算可以elem401(2e11, 0.3, 0.01)
直接返回对应的刚度矩阵
elem401.m
% 四节点单元刚度矩阵函数 E杨氏模量, nu 泊松比, h板厚度
function Ke = elem401(E,nu,h)
syms x
syms y
M = [
[1, x, y, x*y,0,0,0,0]
[0,0,0,0, 1, x, y, x*y]
];
A = [
[1 -1 -1 1 0 0 0 0]
[0 0 0 0 1 -1 -1 1]
[1 1 -1 -1 0 0 0 0]
[0 0 0 0 1 1 -1 -1]
[1 1 1 1 0 0 0 0]
[0 0 0 0 1 1 1 1]
[1 -1 1 -1 0 0 0 0]
[0 0 0 0 1 -1 1 -1]
];
N = M/A; % M * inv(A) , 为形状函数矩阵
% 对x, y等参变量求导并得到对应的应变矩阵
N_x = diff(N,x);
N_y = diff(N,y);
B = [N_x(1,:)
N_y(2,:)
N_y(1,:) + N_x(2,:)];
clear N_x N_y
D = E * h/(1-nu^2) * [
[1 nu 0]
[nu 1 0]
[0 0 (1- nu)/2]
];
k_e = transpose(B) * D * B;
Ke = int(int(k_e,x,[-1, 1]), y, [-1,1]); % 积分求解对应的四节点刚度矩阵
end
另外需要注意的是,单元刚度矩阵一定是奇异矩阵,即
d
e
t
(
K
e
)
=
0
det(K_e) = 0
det(Ke)=0, 如果需要显示矩阵,使用vpa(K_e)
即可,但是注意这个会损失精度导致求解
d
e
t
det
det不为0
下面是八节点单元的matlab代码
调用方法是K_e = Elem801(E,nu,h)
节点编号顺序为
源程序Elem801.m
function Ke = Elem801(E, nu,h)
syms x y
% 对于平面八节点单元,共有
M = [
[1 x y x^2 x*y y^2 x^2*y x*y^2 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 1 x y x^2 x*y y^2 x^2*y x*y^2]
];
% 定义八节点单元的x_e和y_e
x_e = [-1 0 1 1 1 0 -1 -1];
y_e = [-1 -1 -1 0 1 1 1 0];
% 创建初始矩阵A
A = zeros(16,16);
for i = 1:8
A(2*i-1:1:2 *i,:) = subs(M,[x,y],[x_e(i),y_e(i)]);
end
% 计算形状函数矩阵
N = M/A;
% 使用求导方法求解对应的位移函数
N_x = diff(N,x);
N_y = diff(N,y);
B = [N_x(1,:)
N_y(2,:)
N_y(1,:) + N_x(2,:)
];
clear N_x N_y
% 刚度矩阵D
D = E/(1-nu^2) * [
[1 nu 0]
[nu 1 0]
[0 0 (1- nu)/2]
];
ke_pre = transpose(B)*D*B; % 注意不能使用B'
Ke = h * int(int(ke_pre,x,[-1, 1]), y, [-1,1]); % 积分得到刚度矩阵
% 注意不能使用vpa函数, vpa函数会产生截断误差
clear ke_pre
end
% Test Code: Elem801(2e5,0.3, 10)