基于MATLAB梁模型振型动画程序设计

基于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 = [0 0 40 40 80 80;
0 40 0 40 0 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([2 * nodes - 1; 2 * nodes], 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 = [ones(1, 2); localCoords’];
lengthOfElement = abs(det§);
C = inv§;
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 = [1 2];
doffsDirichlet = [Dirichlet * 2 - 1, Dirichlet * 2];
K(doffsDirichlet, ? = 0;
M(doffsDirichlet, ? = 0;
M(doffsDirichlet, doffsDirichlet) = eye(numel(doffsDirichlet));

% 求解固有频率与模态振型
[a, lambda] = 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).