梁模型振型动画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(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 = [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/matlabc ... introduction-to-fem