MATLAB Answers

0

Spectral analysis of 1D elastic wave

Asked by Mohammadfarid ghasemi on 11 Sep 2019
Latest activity Commented on by Mohammadfarid ghasemi on 12 Sep 2019
As indicated in the: "wave propagation in structures: an FFT-based spectral analysis methodology", by James F. Doyle:
The spectral solution of 1D- elastic wave equation is as follow:
where is the wavenumber. C and D are the undetermined amplitudes at each frequency. Let the end of the bar at x = 0 be subjected to a force history F(t), that is, . E and A are the Elastic moduli and cross-sectional area respectively. The final solution is the inverse Fourier transform of the following expression:
is the Fourier transform of the applied force F(t).
The numerical example for the above problem is provided as follow:
Rod:
diameter=1 inch
density=0.00247 lb/ci
E=10.6e6 lb/si
Pulse, F(t):
0.000000 0
0.001000 0
0.001100 1000
0.001300 0
0.001500 0
(sec) (N)
I wrote the following code in MatLab:
clear all
close all
clc
d=1.0; %inch
A=pi/4*d^2;
rho=0.000247; %lb/inch3
E=10.6e6; %psi
%transform parameters:
n=2^15;
dt=5e-6;
fs=1/dt;
time_fcn = (0:n-1)/fs;
frequency = (0:n-1)*(fs/n);
omega=2*pi*frequency;
F=zeros(1,numel(time_fcn));
nn=find(time_fcn>=0.0011 & time_fcn<=0.0013);
F(nn)=-5e6*(time_fcn(nn)-0.0013);
plot(time_fcn,F)
Fn=fft(F);
plot(omega,Fn)
k=omega*sqrt(0.000247/10.6e6);
A=-Fn(2:numel(omega))./(1i*k(2:numel(omega))*E*A);
x=0;
G(2:numel(omega))= A.*exp(-1i*k(2:numel(omega))*x);
G(1)=simpsons(F,0,max(time_fcn),numel(time_fcn));
U= ifft(G);
plot(time_fcn*1000,U)
The result must be as follow:
However, I cannot get the same result as indicated in the abovementioned book. Can anybody tell me where is my mistake?
Thank you all,
Regards.

  2 Comments

darova
on 11 Sep 2019
Here ie written F = 1000 at Pulse = 0.0011
% Pulse, F(t):
% 0.000000 0
% 0.001000 0
% 0.001100 1000
% 0.001300 0
% 0.001500 0
% (sec) (N)
nn = find(time_fcn>=0.0011 & time_fcn<=0.0013);
F(nn)=-5e6*(time_fcn(nn)-0.0013); % 5 000 000 ?
Why don't use constants you defined?
% k = omega*sqrt(rho/E);
k = omega*sqrt(0.000247/10.6e6);
Use additional variables for clearer coding
ii = 2:numel(omega);
% overwrite constant A ?
% A=-Fn(2:numel(omega))./(1i*k(2:numel(omega))*E*A);
FEA = -Fn(ii)./(1i*k(ii)*E*A);
x = 0;
G(ii) = FEA.*exp(-1i*k(ii)*x);
Is it integration?
G(1) = simpsons(F,0,max(time_fcn),numel(time_fcn));
Also use hold on before new plot to add new data
plot(x,y)
hold on
plot(x1,y1)
hold off
Dear Darova,
Thank you for your commnet. I wonder if there is any solution to obtain the same results as in the figure attached.

Sign in to comment.

0 Answers