梁模型振型动画MATLAB程序
%该模型为一端固定梁模型振型动画clear; close all; clc
% 系统参数
E = 1e7;
A = 1.5;
rho = 2.6e-4;%密度
% 节点坐标矩阵第一行为x坐标,对应列第二行为y坐标.既有坐标点(0 0)，(0 40)，(40 0)，(40 40)，(80 0)，(80 40)
p = [00 40 40 80 80;
0 400 400 40];
numberOfNodes = size(p, 2);% 取出系统节点个数即p.
% 对各节点进行建模如第一行表示：连接1 3节点形成一个单元，第二行1 4表示连接1 4节点形成一个单元
t = [1 3;
1 4;
2 4;
3 4;
3 5;
4 5;
4 6;
5 6];
numberOfElements = size(t, 1);% 取出系统总的单元个数即t.
% c
c = A * E;
% Initialization of K and F
K = zeros(2 * numberOfNodes);% 初始刚度阵
M = zeros(2 * numberOfNodes);
F = zeros(2 * numberOfNodes, 1);
% 以下for循环为计算整体刚度阵，质量阵为计算固有频率做准备，经过试验下列for循环可以使用自己写成的有限元整体刚度阵，质量阵来代替。但要注意
%在amp = a(:, mode)' * exp(1i * omega(mode, mode) * timeStep); amp = reshape(amp, 2, numberOfNodes);中输入的numberOfNodes要与自己实际模型一致
%
for e = 1 : numberOfElements
nodes = t(e, :);
dofs = reshape(, 1, 2 * numel(nodes));
nodeCoords = p(:, t(e,:));
n = diff(nodeCoords, 1, 2);
n = n / norm(n);
Q = [n(1) n(2) 0 0;
0 0 n(1) n(2)];
localCoords = Q * nodeCoords(:);
P = ;
lengthOfElement = abs(det(P));
C = inv(P);
diffPhi = C(:, 2);
Ke = Q' * diffPhi * A * E * diffPhi' * lengthOfElement * Q;
% Local shape functions
phi_1 = @(x) C(1,1) + C(1,2) * x;
phi_2 = @(x) C(2,1) + C(2,2) * x;
a = localCoords(1);
b = localCoords(2);
intPhi = [integral(@(x)phi_1(x) .* phi_1(x), a, b) 0 integral(@(x)phi_1(x) .* phi_2(x), a, b) 0;
0 integral(@(x)phi_1(x) .* phi_1(x), a, b) 0 integral(@(x)phi_1(x) .* phi_2(x), a, b);
integral(@(x)phi_1(x) .* phi_2(x), a, b) 0 integral(@(x)phi_2(x) .* phi_2(x), a, b) 0;
0 integral(@(x)phi_1(x) .* phi_2(x), a, b) 0 integral(@(x)phi_2(x) .* phi_2(x), a, b)];
Me = intPhi * rho * A;
K(dofs, dofs) = K(dofs, dofs) + Ke;
M(dofs, dofs) = M(dofs, dofs) + Me;
end
%以上为生成整体刚度阵，质量阵循环程序
% 施加固定边界条件
Dirichlet = ;
doffsDirichlet = ;
K(doffsDirichlet, :) = 0;
M(doffsDirichlet, :) = 0;
M(doffsDirichlet, doffsDirichlet) = eye(numel(doffsDirichlet));
% 求解固有频率与模态振型
= eig(M \ K);% lambda为特征值，a为与特征值对应特征向量，模态振型主要由特征向量来生成
[~, permutation] = sort(diag(lambda));% 对lambda特征值进行从大到小排序，并存入[~, permutation]矩阵中
lambda = lambda(permutation, permutation);
a = a(:, permutation);
% 求解转速
omega = sqrt(lambda);%实模态分析法因此需要对结果开方
% 求固有频率
frequencies = diag(omega / (2 * pi));
sprintf('f = %.1f Hz\n', frequencies)%输出各阶固有频率
% 画图
set(gcf, 'color', 'w')
mode = 5;%确定画第几阶模态振型图
gain = 3.5;
tEnd = 3 * 2 * pi / omega(mode, mode); % 3 periods
for timeStep = 0 : tEnd / 5e1 : tEnd%为下列插值模态动画进行时间步设计
cla
amp = a(:, mode)' * exp(1i * omega(mode, mode) * timeStep);%求个时间步下的模态位移变化量
amp = reshape(amp, 2, numberOfNodes);%为将变化后的模态位移变化量加到原模型上，对结果矩阵进行一定变换
pnew = p + gain * imag(amp);%将模态位移变化量加到原来的模态位移之上
hold on
for e = 1 : numberOfElements
nodes = t(e, :);
plot(p(1, nodes), p(2, nodes), 'b--');
plot(pnew(1, nodes), pnew(2, nodes), 'k-o',...
'MarkerFaceColor', 'k', 'MarkerSize', 10, 'LineWidth', 3);
end
axis equal
axis([-5 85 -5 45])
axis off
drawnow
end
源程序：https://ww2.mathworks.cn/matlabcentral/fileexchange/46384-introduction-to-fem
页:
[1]