[nikla]

There was conducted ten different experiments, with ten different velocities in an air chamber where air at 25 degrees Celsius was blown at a cylinder of 20mm. The measurements were taken by a hot-wire anemometer a certain distance behind the cylinder. The frequency of measurements were 1000Hz and we were given ten different files with one array in them with velocities.

The end goal is too make two graphs, one showing vortex shed freq vs Velocity, and Strouhal- vs Reynolds number as a function of Re, and finding uncertainties considering dU=0,1. I do know the formulas for the two last variables:

St = f*D/U
Re = rho*U*D/mu

But I am not sure how to get the shedding frequency. And if you look at the image under, I am also unsure how to deal with the very start of the FFT.


I am not sure if I am supposed to do anything with the start like windowing or what not.


What I have done here is to use FFT and nyquist limits and making PSD out of it. I have the impression that if i take the max value of each of these plots and plot them against the velocity I will have what i requested first. However with the start value being so high for each of them the plot is a straight line.

I am not looking for anyone to do it for me, but rather if anyone could point me to relevant information about my problem. If it was not clear i am doing this in MATLAB.

[FMDenaro]

First of all, the zero-th wavenumber in the spectra is the average energy content, try to plot the spectra in a log-log scale without the zero.
Then, consider that your time-signal is not exactly periodic, are you performing the FFt over the whole interval or are you using windowing technique?

[nikla]

Originally Posted by FMDenaro

First of all, the zero-th wavenumber in the spectra is the average energy content, try to plot the spectra in a log-log scale without the zero.
Then, consider that your time-signal is not exactly periodic, are you performing the FFt over the whole interval or are you using windowing technique?

Thanks for the reply!


I am not currently using a windowing technique, as I am unsure if its appropriate in this case.


When I should remove the zero-th wave number is this before or after the FFT? If it is the later i get this for my plots:





Do these to plots look like something that could be correct?

[FMDenaro]

Remove only from the plot, first of all show the energy spectra resulting from the whole signal in a log-log plot

[nikla]

Originally Posted by FMDenaro

Remove only from the plot, first of all show the energy spectra resulting from the whole signal in a log-log plot

Okey so I plotted this:

loglog(PSDx{i})

From:

%% Number of points
    N = numel(Ucell{i}); % Lenght of signal
    
    %% FFT
    xdft{i} = real(fft(Ucell{i}));
    xdft{i} = xdft{i}(1:N/2+1); % Nyquist
    
    %% PSD - Power Spectral Density
    PSDx{i} = (1/(fs*N)) * abs(xdft{i}).^2;
     PSDx{i}(2:end-1) = 2*PSDx{i}(2:end-1);

I then get this:

[FMDenaro]

what are you exactly representing in x and y axes? If 10^3 Hz is the Nyquist cut-off how can you plot higher frequencies??? What is exactly the time-step of your measurement device?

[nikla]

Originally Posted by FMDenaro

what are you exactly representing in x and y axes? If 10^3 Hz is the Nyquist cut-off how can you plot higher frequencies??? What is exactly the time-step of your measurement device?

Sorry I tought wrong, and plotted wrong.

It now cuts of at the nyquist limit 500hz.

%% Load test data
Ucell   = {load('3d7.lvm') load('5d5.lvm') load('7d1.lvm')...
    load('8d7.lvm') load('11d9.lvm') load('14d4.lvm') load('16d5.lvm')...
    load('18d9.lvm') load('20d8_1.lvm') load('20d8_2.lvm')};
UL      = numel(Ucell);
U       = [3.7 5.5 7.1 8.7 11.9 14.4 16.5 18.9 20.8 20.8]; % m/s

%% Variables and Conditions
D       = 0.02;                 % [m] Diameter of cylinder
T       = 25;                   % [C] Tempratur
rho     = 1.184;                % [kg/m^3] Air density
mu      = 18.37e-6;             % [N*s/m^2] Dynamic (absolute) viscosity
deltaU  = 0.1;                  % [m/s] speed change
fs      = 1000;                 % [Hz] Sampling frequency
ts      = 1/fs;                 % [s] Sample time
Ns      = fs/2;                 % [Hz] Nyquist sampling frequency

%% Calculating
% Preallocation

for i=1:UL
    %% Number of points
    N = numel(Ucell{i}); % Lenght of signal
%     t{i} = (0:N-1)*ts; % Time vector
    
    %% FFT
    xdft{i} = real(fft(Ucell{i}));
    xdft{i} = xdft{i}(1:N/2+1); % Nyquist
    
    %% PSD - Power Spectral Density
    PSDx{i} = (1/(fs*N)) * abs(xdft{i}).^2;
    PSDx{i}(2:end-1) = 2*PSDx{i}(2:end-1);
    freq{i} = 0:fs/N:Ns;
    
    %% Max value
    [M(i),I(i)] = max((PSDx{i}));
    I(i) = (I(i)*Ns)/N; % Getting the frequency and not indic place

    %% Strouhal number and Reynolds number
    St(i) = (M(i)*D)/U(i); % Strouhal number
    Re(i) = (rho*U(i)*D)/mu; % Reynolds number
end


%% PSD
figure(4)
for i=1:UL
    subplot(numRow,plotsPerRow,i)
    loglog(freq{i},PSDx{i})
    title(sprintf('U=%0.5g [m/s]', U(i)))
     xlabel('Frequency (Hz)')
     ylabel('PSD')
    grid on
end
sgtitle('Loglog')
linkaxes % Same axis range for all

[FMDenaro]

But you must not cut the plot, your whole signal contains frequencies only up to the Nyquist one and that means you perform the FFT on the whole signal and then plot the squared modulus of all Fourier couefficients along each frequency.

And the Nyquist frequency is pi/dt, what is yout dt?