Spectral analysis of 1D elastic wave

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]

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]

Dear Darova,

Thank you for your commnet. I wonder if there is any solution to obtain the same results as in the figure attached.

[Foil]

Hi folk，

Have you solved your problem?

I'm looking for a code to simulate the elastic wave propagation in an rod or plate or frame,...!

Would you please help me to find such a thing.

Regards.