MATLAB Answers

Spectral analysis of 1D elastic wave

10 views (last 30 days)
Mohammadfarid ghasemi
Mohammadfarid ghasemi on 11 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
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
Mohammadfarid ghasemi
Mohammadfarid ghasemi on 12 Sep 2019
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.

Answers (0)

Sign in to answer this question.


Translated by