mechanics of vibration

Image_input_New.m

clc clear all A = imread('lute1.gif'); RangeX = size(A, 2); RangeY = size(A, 1); A = double(A)/double(max(A(:))); B = (1-(A == 1)); [UnitY, UnitX] = find(B); Border_Ind = boundary(UnitX, UnitY, 1); Border_X = UnitX(Border_Ind); Border_Y = UnitY(Border_Ind); Body_length = 0.7; % length in meter Unit_length = Body_length / RangeX; a = Unit_length % unit length of lump mass, in meter E = 11.03e9; % Elastic modulus in Pa Density = 0.45e3; % Density of Spruce wood (for top board), in kg/m3 thickness = 0.005; % Thickness in meter Alpha = 0.0611; % Square board bending constant Unit_Mass_board = Density * thickness * a^2; Unit_Mass_side = Unit_Mass_board * 30; K_equi = 3.0023e-4 ; % K_equi = K * Alpha * a^2 / E / thickness^3 K_board = K_equi / Alpha / a^2 * E * thickness^3; K_side = K_board * 1000; MAP = zeros (max(UnitY), max(UnitX)); N = length(UnitY); M = Unit_Mass_board * ones(1, N); Link = zeros(N, N); Link_Self = zeros(1, N); K_link = zeros(N, N); KK = zeros(N, N); XX = zeros(N, N); w = zeros(N, N); Diff_Y = zeros(N, N); Diff_X = zeros(N, N); X = UnitX; % * Unit_length; Y = UnitY; % * Unit_length; for j = 1:N for k = 1:N Diff_X(j, k) = abs(UnitX(j)-UnitX(k)); Diff_Y(j, k) = abs(UnitY(j)-UnitY(k)); end end Link = ((Diff_X + Diff_Y) == 1); Link_in_X = Link .* Diff_X; Link_in_Y = Link .* Diff_Y; for j = 1:N K_link (j,:) = Link_in_X(j,:) * K_board + Link_in_Y(j,:) * K_board; Link_Self(j) = sum(Link(j, :)); end for j = 1:N KK (j,:) = - K_link (j,:); KK (j,j) = - sum(KK (j,:)); end Edge_Ind = find(Link_Self < 4); Hole_Ind = []; % elements that are next to a hole Side_Ind = setdiff(Edge_Ind, Hole_Ind); % Exclude hole elements from edge elements and call the rest side elements side elements B(Hole_Ind) = 2; for j = 1:length(Side_Ind) k = Side_Ind(j); KK(k,k) = KK (k,k) + K_side; end M(Side_Ind) = Unit_Mass_side; MM = diag(M); [XX, w_sqrt] = eig(KK, MM); Nat_Freq = ones(1, N) * (w_sqrt.^0.5); save Vibration_without_Hole MM KK X Y Nat_Freq XX MAP Side_Ind Hole_Ind Border_X Border_Y UnitX UnitY for j = 1:N MAP(UnitY(j),UnitX(j)) = XX(j,5); end loops = 25; Frame(loops) = struct('cdata',[],'colormap',[]); PP = surf (MAP,'FaceAlpha',0.5); zmax = max(max(abs(MAP))); axis([min(X) max(X) min(Y) max(Y) -zmax zmax]); colormap default hold on plot3 (Border_X, Border_Y, zeros(length(Border_X)), '-', 'LineWidth', 1.5, 'color', 'r') hold off ax = gca; af = gcf; ax.DataAspectRatio = [1,1, zmax/max(X)* 10]; af.Position = [100, 100, 1000, 600]; for j = 1:loops t = 2*pi/loops * (j-1); MAP1= MAP.*cos(t); set(PP, 'ZData',MAP1); drawnow; Frame(j)=getframe; end set(gcf, 'Position', [100, 100, 800, 500]);movie(gcf,Frame, 50, 20)