Spectral analysis of 1D elastic wave

Mohammadfarid ghasemi (view profile)

on 11 Sep 2019
Latest activity Commented on by Mohammadfarid ghasemi

Mohammadfarid ghasemi (view profile)

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.

darova

darova (view profile)

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