Dispersion Curves Flexural And Longitudinal Waves

% Analysis of Wave Propagation in Structures and Solids
%
% Exercise 4 (FEM and WFE)

clear all
close all
clc

%% Definition of the bar

Le=0.006; % length of the beamelement [m]
rho=7800; % Density [kg/m^3]
h=0.0005; % heigth [m]
E=2e11; % Young modulus [Pa]
b=0.0005; % width [m^2]
nu=0.329;% Poissons ratio


I=(b*h^3)/12; % Second moment of area of cross-section about the neutral axis
omega=0:5:2e5; % Defining omega from 0 to 5e5 in steps of 5.%% The problem will be solved twice for longitudinal and flexural waves.
% Subindex l denotes longitudinal waves whereas f denotes flexural waves.

%% Analytical solutions (wavenumber for thedifferent types of waves)

kl=omega.*(rho/E)^(1/2); % longitudinal wavenumber
kf =((omega.^2.*rho*h*b)/(E*I)).^(1/4); % flexural wavenumber

%% WFE solution

elem=1; %Always1 element for WFE
a=Le/(2*elem); %Half of the element

%% Calculation of the dispersion curves for the longitudinal waves

%Mass and stiffness matrices for longitudinal waves (givenin previous exercises)

Kl= (E*h*b)/(Le)*[1, -1; -1, 1];

Ml= (rho*h*b*Le)/(6)*[2, 1; 1, 2];

lambdal=exp(-i*kl*Le);

Il=eye(size(Kl)/2); % Identity matrix for longitudinalwaves' analysis (size half of the Kl)

WFEomegal=zeros(size(kl));

for ii=1:length(lambdal)
LambdaL=[Il 1/lambdal(ii)*Il];
LambdaR=[Il, lambdal(ii)*Il].';
Kbarl=LambdaL*Kl*LambdaR;Mbarl=LambdaL*Ml*LambdaR;
sol=eig(Mbarl\Kbarl);
sol=sort(sol);
WFEomegal(ii)=real(sqrt(sol(1)));
end

%% Calculation the dispersion curves for the flexural waves

%Mass andstiffness matrices for flexural waves (given in the lecture notes)

Kf = ((E*I)/(2*a^3))*[3, 3*a, -3, 3*a;
3*a, 4*a^2, -3*a, 2*a^2;
-3, -3*a, 3, -3*a;
3*a, 2*a^2, -3*a, 4*a^2];

Mf =...