Simulation numérique de l’instabilité vibratoire à nombre de Reynolds élevé dans un tandem de cylindres et un faisceau de tubes
Christophe HEUDES & Jingyan DUO
Supervisor : Marianna Braza and EDF
Numerical simulation of vibrational instability with high Reynolds
number on cylinders in tandem and a bundle of tubes
Introduction
During our project, we have worked on fluid-structure interactions on two cylinders in tandem with numerical simulation. Structures in contact with fluid flow, whether natural or man-made, are subject to flow-induced forces and flow-induced vibration. Under certain conditions, the vibration may be self-excited. These instabilities cause damage in a short or the long term in nuclear power plants for example. It is why the BARESAFE project was created (it regroups EDF, AREVA and CEA) to understand these instabilities. Firms have to replace cylinders every six months due to the instabilities ; it is clearly important to understand and fix this problem. Therefore our main interest is twofold because these cylinders can be found in offshore installations or cables.
Our aim is to build a reduced model and a post-process with a big data based on an experience performed in a wind tunnel to decrease calculation times. Indeed, if calculation times are reduced, we can test more complex configurations and reach as quickly a resolution as possible.
Présentation du problème
To study the problem, we will consider two cylinders in tandem which is the simplest form. The configuration is the same as what is been studied by NASA in their experiments so we can easily compare our results with the experimental results. The first configuration was used to study noise of a aircraft because it modeled the supports of the aircrafts’ landing gear. We will use the same Reynolds 166 000 as a first step.
Figure 1.1 – Cylinders in tandem
Fluid-structure interaction
Consider that we have a characteristic length D with a natural frequency fs in a flow which velocity is U. Dimensional analysis enables to show a dimensionless number u* which is the reduced velocity.
This parameter enables us to describe the fluid-structure interaction phenomena, it compares the characteristic time ofthe structure and the characteristic time of advection. By plotting amplitude of oscillation as function of the reduced velocity, we can consider three zones which correspond to different mechanisms :
Figure 1.2 – Amplitude of oscillation as function of reduced velocity
— Zone 1 corresponds to Turbulence Induced Vibrations (TIV). The amplitude is limited and increases steadily with the reduced velocity.
— Zone 2 is the zone of Vortex Induced Vibrations (VIV). It looks like a resonance where natural modes are excited by the wake. Amplitude increases but the system remains stable.
— Zone 3 corresponds to Movement Induced Vibrations (MIV), the system may be unstable and cause the break of the structure. 2
Modélisation numérique
NSMB Code
In 1989, Navier Stokes Multi-Block (NSMB). This code has been created by european research laboratories in aerodynamicsuch as École Polytechnique Fédérale de Lausanne (EPFLausanne), KTH, CERFACS, IMFT, EADS, SAAB Military Aircraft, CFS Engineering, ETH-Zurich, IMFS de Strasbourg. Its was designed for parallel computing by solving blocks in parallel on different processors. NSMB enables to solve Navier-Stokes equations on a mesh divided in several blocks. Thanks to this approach, we can mesh the structure with a complex geometry
Simulation
Figure 2.1 – Problem scheme
During the first week, we worked on a 2D model SST k-omega. As a first step, we worked on a simulation where the two cylinders were immobile to study results on position cylinder and force coefficients.
Turbulence models
Averaged equations for turbulence model are :
Boussinesq approximation is used to model Reynolds stress tensor :
It is a model with two equations on k and w :
For DDES kw-sst the model length scale is calculated as :
Model OES (Organised Eddy Simulation) is used when turbulence is not in equilibrium and enables to see Von Kármán effects. This model is based on the k epsilon model. Equations for k epsilon model are :
2D simulation
Parameters chosen for the 2D simulation are :
— Mach number 0.1285
— Reynolds 166 000
— Gas constant 43.25
— Timestep 0.00845
— CFL 1
— residue variable : density
— Tolerance 1e − 4
— Time scheme : implicit
— space scheme : centered
Eventually, we used a mesh where the second cylinder is free, enabling it to move axially upwards. This configuration represents the movement of cylinders in nuclear power plants in the coolant circuit. The following animation is what we have obtained with Tecplot, during one period ( 1000 snapshots ), with u* = 4 which is a supercritical regime.
3D simulation
For the 3D simulation, two turbulent models were given to us : OES model and DDES k omega-SST model. 
Figure 2.2 – 3D Visualization with kw SST model
We have tested these two models to study the differences, the aim is to perform the 2D study on these 3D cases and follow the Kelvin-helmholtz structures.
Figure 2.3 – Slice on q criterion on the 3D model
We could represent the pressure coefficient around the second cylinder by using :
Then, we can plot the pressure coefficient around the second cylinder compared with NASA results :
Figure 2.4 – Pressure coefficient around the second cylinder compared with NASA results
Analysis of the results
Landau global oscillator model
In this part, we will study the Landau global oscillator model with a Re =132812.
Once we have the data of the vibration position, we will use the Landau global oscillator model to analysis the physical parameters. The amplification envelopes of the signals have been treated by the Landau global oscillator model (Landau, 1944), Provansal and al (Provansal and al., 1987). Here is the Stuart-Landau equation:
\begin{equation}
\frac{dA}{dt}=\sigma A-\frac{1}{2} l |A|^2 A
\end{equation}
Where A is the complex amplitude of the transverse position in
the flow axis, which is the most significant quantity of the symmetry (varies on time especially for the transition regime); $\sigma$ is the relative growth rate; $l$ is the Landau constant.
Then we have the Landau equation :
\begin{equation}
\frac{d|A|}{dt}=\sigma_r |A|-\frac{1}{2} l_r |A|^2 A
\end{equation}
where $\sigma_r$ is the amplification rate in the linear regime and $l_r$ is the Landau constant in the non-linear regime.
Critical reduced velocity $u^*_{critical}$
In this part, we are going to calculate the critical value of the reduced velocity. At first, we will have a diagram of $log(A)$ as a function of real time, A is the amplitude of the vibration of the cylinder. In our study, the position of the cylinder is dimensionless but the calculate time is real time.
Figure 3.1 –Different regime of the a signal, example $u^*=4$
In the linear regime, the predominant term is $\sigma_r$; in the non-linear regime, the predominant term is the $l_r$. In other words, in the following study, in order to calculate the $\sigma_r$ we will analyze the signal in the linear regime; on the contrary, the signal in the non-linear regime is more interesting for us to study the Laudau constant $l_r$.
Once we have the signal data, we should at first identify the linear regime and then draw the diagram of $log(A)$ as a function of the period. We will do the same process to have the function of the fitting line for each case of reduced velocity, especially the slope of these fitting lines, amplification rate in the linear regime $\sigma_r$. In order to have the value of the critical reduced velocity, we should draw the $\sigma_r$ as a function of $u^*$, the intersection between the fitting line and $\sigma_r=0$ indicates the value of the critical reduced velocity $u^*$.
Figure 3.2 – Critical reduced velocity $u^*$
With this study, we have obtained $u^*_{critical}=3.55$ when $\sigma_r=0$
Landau constant $l_r$
In order to calculate the Landau constant, we should at first identify the transition regime and the the saturated regime, where the non-linear term $\frac{1}{2} l_r |A|^2 A$ dominates in the Landau equation. With the diagram of $A$ as function of $T A^2$, we can have the slope of the fitting line of these points.
Figure 3.3 – Landau constant $l_r$ for $u^*=3.3$
For the case $u^*=3.3$ we have obtained that the Landau constant $l_r=-5.2<0$, which indicates a sub-critical regime, since here we have $u^*=3.3<u^*_{critical}$ ($u^*_{critical}=3.55$).
Spectral analysis
In this part we will analysis the Strouhal number and the kelvin Helmholtz instability by the frequency. We have defined 8 monitor points in order to calculate the evolution of the physical parameters on the defined points ( especially the pressure).
Here we have chosen a horizontal line, since the wake moves periodically and alternatively, so the horizontal line is the average position of the wake.
Figure 3.4 – Position of the monitor points
In the TABLE 3.1. are the coordinates of the eight monitor points. The block and the coordinate (k, i and j) are defined in the mesh reference of the NSMB.
| monitor point | block | k | i | j |
| 1 | 13 | 1 | 83 | 12 |
| 2 | 13 | 1 | 64 | 26 |
| 3 | 13 | 1 | 55 | 30 |
| 4 | 13 | 1 | 45 | 36 |
| 5 | 13 | 1 | 39 | 41 |
| 6 | 5 | 1 | 6 | 22 |
| 7 | 5 | 1 | 11 | 24 |
| 8 | 5 | 1 | 20 | 28 |
Table 3.1 – Coordinates of the monitor points
Then we will use the FFT approach to observe the variation of the PSD as a function of the St ( dimensionless value of the frequency ).
FFT approach estimate
As there are several methods of FFT approach method, we will at first identify which method works better for our study. Here we will talk about two methods, Periodogram PSD estimate and Welch PSD estimate. With several calculations, we have found that, periodogram PSD estimate works well for the low frequency signal; Welch PSD estimate has a good estimation for the high frequency signal. So in our study, we will fix a critical valueof the frequency $f_{cri}$ ( $St_{cri}$ ), if $f<f_{cri}$ ( $St<St_{cri}$ ), we will use the periodogram PSD estimate; if $f>f_{cri}$ ( $St>St_{cri}$ ), we will use the Welch PSD estimate. In the following charts, we can clearly identify the Von Kármán vortex shedding which corresponds to a maximum value of the PSD ( the first peak in the diagram ). Based on the study of M.Elhimer (2014), the experimental Kelvin Helhomltz instability frequency value is around 1800 Hz, on which we will have a discuss in detail in the following part.
Figure 3.5 – Variation of PSD as a function of St of monitor ponit 1 on the first cylinder on 2D
In the chart, we find that the Strouhal number equals to 0.23 which matches well our previous calculate ( $St_K=0.23$ ).
Figure 3.6 – Variation of PSD as a function of frequency of monitor point 1 on the first cylinder on 2D
In this chart we find that the Von Kármán vortex frequency $f_K=174.22 Hz$ and we have as well marked the experimental Kelvin Helhomltz instability frequency value which is $f_{SL}=1800 Hz$.
In these two pictures, the black lines represent the final PSD estimate : equals to green line when $f<f_{cri}$ for the low frequency signal; equals to red lines when $f>f_{cri}$ for the high frequency signal. In both pictures, the results of the Von Kármán shedding matches well with our previous study, we have marked the experimental values for the Kelvin Helhomltz instability with $f=1800 Hz$ and $St=2.3$. We have also performed this study for the other monitor points and also for the 2nd cylinder, the results are almost the same.
Time-frequency-transformations
CWT-Continuous wavelet transform
In mathematics, a continuous wavelet transform (CWT) is used to divide a continuous-time function into wavelets. Unlike Fourier transform, the continuous wavelet transform possesses the ability to construct a time-frequency representation of a signal that offers very good time and frequency localization.
Figure 3.7 – Wavelet transform of the 4th monitor point on the 1st cylinder on 2D
In this picture, in order to observe the Kármán vortex frequency, which corresponds to a $St=0.23$ ( 514 points in a period ), we have use a defined scales from 450 to 600. We have found exactly the same value of the $St_K=0.23$.
AR-Autoregressive power spectral density estimate
Compared to the FFT approach, AR estimate has a better stability for short segments of signal, a better spectral resolution and a better resolution as function of time. We will present the results on the 4th monitor point in the following part.
Figure 3.8 – AR transform of the 4th monitor point on the 1st cylinder on 2D
Here we have observed the $St_K=0.23$ which corresponds well to our previous study. But in this case, we can't observe the phenomena of shear layer vortex. So in the following part, we will at first filter the low frequency signal ( for example, $Fc=2.2$ for the strouhal number ), then we excuter a AR transform.
Figure 3.9 – High-pass filtered signal of the variation pressure
Here we define the $Fc=2.20$ in order to observe the Strouhal number corresponding to the Kelvin Helhomltz instability.
Figure 3.10 – AR transform of the 4th monitor point on the 1st cylinder on 2D with a high-pass filter Fc=2.2
With a high-pass filter, we can observe the Kelvin Helhomltz instability phenomena which corresponds to a $St=2.3$.
The method of Proper Orthogonal Decomposition
The Proper Orthogonal Decomposition (POD) is a post-processing technique. It is a multi-variate statistical method that aims at obtaining a compact representation of the data. This method may serve two purposes, namely order reduction by projecting high-dimensional data into a lower-dimensional space and feature extraction by revealing relevant, but unexpected, structure hidden in the data.
It takes a given set of data and extracts basis functions, that contain as much “energy” as possible. By taking first modes of POD, we can modelize almost all the problem because they have all the energy.
Let $ q_{k}(x), k = 1 ... N_t $ be a set of observations (called snapshots) at point x of the domain that could be obtained by a PIV experimental measurements or by numerical simulation. The goal of POD is to find functions $\Phi$ such that :
$$ \frac{<\mid(q,\Phi)\mid^{2}>}{\mid\mid\Phi\mid\mid^{2}} $$
is maximized. $<.>$ is an average, $\mid\mid.\mid\mid$ the induced norm. Solving the optimization problem leads to an eigenvalue problem, where the functions $\Phi$ are the eigenfunctions.
In our problem, the resolution gives us the eigenfunctions and by taking first modes, we can modelize the velocity :
Figure 3.11 - First mode for velocity corresponds to the mean (numerical results)
Figure 3.12 – POD first modes
Figure 3.13 – POD reconstruction
Conclusion
In general, the numerical simulation is a powerful tool for research as well as for industial engineering. This study allows us to explore a new calculation code NSMB and the parallel computing on supercomputers. We have also studied the Landau theory, different PSD estimate method as well as the POD method in order to calculate the Strouhal number, to observe the phomenema of Von Kármán vortex as well as the Kelvin Helmholtz instability. Future research would be the dimensioning / design of the cylinders bundle in the nuclear reactor based on high fidelity modeling and reduced order/ POD based modeling (BARESAFE ANR project).
Acknowledgements
We express our warm thanks to Mrs Braza, Mr Szubert and Mr Asproulias for their support and guidance during the project. We also thank Mr Deloze and Mr Hoarau for their complete our simulations.
References
[1] S. Bourdet : Analyse physique d’écoulements compressibles instationnaires autour de structures portantes dans le contexte d’interaction fluide-structure. March 2015.
[2] S. BOURDET and al. : Direct Numerical Simulation of the three-dimensional transition to turbulence in the transonic flow around a wing. J. Flow, Turbulence and Combustion, 2003, Vol. 78.
[3] V. Shinde and al. : Numerical simulation of the fluid–structure interaction in a tube array under cross flow at moderate and high Reynolds number, February, 2014.
[4] M. Elhimer and al. : Coherent and turbulent processes in the bistable regime around a tandem of cylinders including reattached flow dynamics by means of highspeed PIV, aInstitut de Mecanique des Fluides de Toulouse, IMFT-UMR CNRS 5502, France.
[5] F. GROSSI : Physics and modeling of unsteady shock wave/boundary layer interactions over transonic airfoils by numerical simulation, Institut National Polytechnique de Toulouse (INP Toulouse), 2014. 20
APPENDIX
Here is our matlab code in order to study data from NSMB (FFT, Welch method...)
%---------------------------------------------------------------------------------------------------------------------------------------------------------
%----------------------------------------------------------------FFT Approach--------------------------------------------------------------------
%---------------------------------------------------------------------------------------------------------------------------------------------------------
% FFT approach for singal
% Author: Jingyan DUO-Christophe HEUDES, ENSEEIHT
% Creation: mars 2015
% Last modification: J. DUO, C. HEUDES
clear all
close all
%PARAMETERS
D=0.05715; %diameter
Uinf=44; %velocity upstream
%fft on pressure data which is the 11th column
display('Please select the file')
monitor=uigetfile
monitor=importdata(monitor);
n=length(monitor(:,11));
%---------------------------------------------------------------periodogram PSD estimate---------------------------------------------------
[psdper,Stt]=periodogram(monitor(:,11)-mean(monitor(:,11)),hann(n),2^(nextpow2(n)+6),1./(monitor(2,3)-monitor(1,3)));
frequence=Stt*Uinf/D;
[y x]=max(psdper)
St_mood1=Stt(x);
f_mood1=St_mood1*Uinf/D;
figure(1)
loglog(Stt,psdper,'k');
hold on;
line([St_mood1 St_mood1],[10e-20 1]);
hold off;
axis([1E-4 1E2 1e-20 1E3]);
title('periodogram PSD estimate as function of St on the monitor point 1 of the 1st cyliner','fontsize',14)
xlabel('Strouhal number','fontsize',12)
ylabel('PSD (Pa?/s)','fontsize',12)
legend('PSD as function of St with periodogram estimate','Strouhal number');
figure(2)
loglog(frequence,psdper,'k');
hold on;
line([f_mood1 f_mood1],[10e-20 1]);
line([1800 1800],[10e-20 1]);
hold off;
axis([1E-2 1E5 1e-20 1E3]);
legend('PSD as function of frequency','Kármán vortex frequency','shear layer vortex frequency')
title('periodogram PSD estimate as function of frequency on the monitor point 1 of the 1st cyliner','fontsize',14)
xlabel('frequency (Hz)','fontsize',12)
ylabel('PSD (Pa?/s)','fontsize',12)
%----------------------------------------------------------------Welch PSD estimate----------------------------------------------------------
[psdwelch1,St] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(2000),1600,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));
[psdwelch2,St] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(1500),1200,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));
[psdwelch3,St] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(1000),800,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));
% %length of 2000 samples with 1700 overlapped samples etc. since here the fs=118.34
% zeropadding with 2^(nextpow2(n)+6)
frequence=St*Uinf/D;
[y x]=max(psdwelch1)
St_mood1=St(x);
f_mood1=St_mood1*Uinf/D;
figure(3)
loglog(St,psdwelch1,'g');
hold on;
loglog(St,psdwelch2,'k');
loglog(St,psdwelch3,'r');
line([St_mood1 St_mood1],[10e-20 1]);
hold off;
xlabel('Strouhal')
ylabel('PSD on the monitor ponit 1')
legend('PSD as function of St with nfft=2000','PSD as function of St with nfft=1500','PSD as function of St with nfft=1000','Strouhal number St=0.23');
title('Welch PSD estimate as function of St on the monitor point 1 of the 1st cyliner','fontsize',14');
figure(4)
loglog(frequence,psdwelch1,'g');
hold on;
loglog(frequence,psdwelch2,'k');
loglog(frequence,psdwelch3,'r');
line([f_mood1 f_mood1],[10e-20 1]);
line([1800 1800],[10e-20 1]);
hold off;
legend('PSD as function of frequency with nfft=2000','PSD as function of frequency with nfft=1500','PSD as function of frequency with nfft=1000','Kármán vortex frequency f=174.13 Hz','shear layer vortex frequency f=1800Hz')
title('Welch PSD estimate as function of St on the monitor point 1 of the 1st cyliner','fontsize',14)
xlabel('frequency (Hz)','fontsize',12)
ylabel('PSD (Pa?/s)','fontsize',12)
%--------------------------------------------------------------------comparaision-----------------------------------------------------------------
%point number 1 of 1st and 2nd cylinder
display('Please select the file of the 1st ponit of 1st cylinder')
monitor1=uigetfile
monitor1=importdata(monitor1);
display('Please select the file of the 1st point of 2nd cylinder')
monitor2=uigetfile
monitor2=importdata(monitor2);
display('Please select the file of the 8th point of the 1st cylinder')
monitor3=uigetfile
monitor3=importdata(monitor3);
display('Please select the file of the 8th point of the 2nd cylinder')
monitor4=uigetfile
monitor4=importdata(monitor4);
figure(5)
plot(monitor1(:,5));
hold on;
plot(monitor2(:,5));
plot(monitor3(:,5));
plot(monitor4(:,5));
hold off;
legend('1st point of 1st cylinder','1st point of 2nd cylinder','8th point of 1st cylinder','8th point of 2nd cylinder');
title('horizontal speed of the flow u/Uinf');
xlabel('step');
ylabel('horizontal speed/Uinf');
%----------------------------------------------------------combinaision des deux methodes----------------------------------------------
clear all
close all
%PARAMETERS
D=0.05715; %diameter
Uinf=44; %velocity upstream
%importer des donness
%fft on pressure data which is the 11th column
display('Please select the file')
monitor=uigetfile
monitor=importdata(monitor);
n=length(monitor(:,11));
%pour bas frequence(f<7699), on utilise le methode de periodogram estimate
%pour haut frequence(f>7699), on utilise le methode de welch PSD estimate
N=10359;%pour amont, point 40
n1=88606;%logeur de bas frequence
n2=44304;%logeur de haut frequence
%diagramme de PSD en fonction de St
[psdper,St1]=periodogram(monitor(:,11)-mean(monitor(:,11)),hann(n),2^(nextpow2(n)+6),1./(monitor(2,3)-monitor(1,3)));
[y x]=max(psdper)
St_mood1=St1(x);
[psdwelch,St2] = pwelch(monitor(:,11)-mean(monitor(:,11)),hann(1000),800,2^(nextpow2(n)+5),1./(monitor(2,3)-monitor(1,3)));
figure(6)
loglog(St1,psdper,'g.');
hold on;
loglog(St2,psdwelch,'r.');
psd1=psdper(1:n1,:);
psd2=psdwelch(n2:length(St2),:);
loglog(St1(1:n1,:),psd1,'k');
loglog(St2(n2:length(St2),:),psd2,'k');
line([10 10],[10e-20 1e4]);
line([St_mood1 St_mood1],[10e-20 1e4]);
line([2.34 2.34],[10e-20 1e4]);
hold off;
title('PSD estimate as function of St on the monitor point 1 of the 1st cyliner, with St_{cri}=10');
legend('Periodogram PSD estimate','Welch PSD estimate','Periodogram PSD estimate for low frquency','Welch PSD estimate for high frequency');
xlabel('Strouhal nomber St');
ylabel('PSD estimate (Pa^2/s)');
%diagramme de PSD en fonction de f
figure(7)
loglog(St1*Uinf/D,psdper,'g.');
hold on;
loglog(St2*Uinf/D,psdwelch,'r.');
psd1=psdper(1:n1,:);
psd2=psdwelch(n2:length(St2),:);
loglog(St1(1:n1,:)*Uinf/D,psd1,'k');
loglog(St2(n2:length(St2))*Uinf/D,psd2,'k');
line([7699 7699],[1e-20 1]);
line([1800 1800],[1e-20 1]);
line([St_mood1*Uinf/D St_mood1*Uinf/D],[1E-20 1]);
hold off;
title('PSD estimate as function of frequency on the monitor point 1 of the 1st cyliner, with f_{cri}=7699 Hz');
legend('Periodogram PSD estimate','Welch PSD estimate','Periodogram PSD estimate for low frquency','Welch PSD estimate for high frequency');
xlabel('frequency (Hz)');
ylabel('PSD estimate (Pa^2/s)');
%---------------------------------------------------------------------------------------------------------------------------------------------------------
%-------------------------------------------------------------------Wavelet transform-----------------------------------------------------------
%---------------------------------------------------------------------------------------------------------------------------------------------------------
% Compute CWT - Simple version
% Author: Damien SZUBERT - IMFT
% Creation: June 2013
% Last modification: D. SZUBERT 23II2015
%
% input :
% - tt: time vector (constant time step),
% - sig: signal (don't forget to remove mean!),
% - wname: wavelet name,
% - scales: scales of CWT (use scal2frq)
% - plotflag: plot type ('surf','imagesc','contourf')
%
% output :
% - cwt_data: coefficient of wavelet tranform,
% - pseudofreq: frequencies of CWT
%
% RMK: scales corresponds to periods in terms of samples (more or less)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [cwt_data,pseudofreq] = monitorCWT_simple(tt,sig,wname,scales,plotflag)
pseudofreq = scal2frq(scales,wname,tt(2)-tt(1)); % Calculate frequency values from scales
cwt_data = cwt(sig,scales,wname); % Compute Continuous Wavelet Transform
%save('cwt_data.mat','cwt_data','tt','pseudofreq','-v7.3')
cwt_data_min = min(min(abs(cwt_data)));
cwt_data_max = max(max(abs(cwt_data)));
disp(['Min / max : ',num2str(cwt_data_min),' / ',num2str(cwt_data_max)])
%% Plot
if strcmp(plotflag,'surf')
surf(tt,pseudofreq,abs(cwt_data),'edgecolor','none');view(2)
elseif strcmp(plotflag,'imagesc')
imagesc(tt,pseudofreq,abs(cwt_data))
warning('Pay attention to the frequency scale! (log/linear axis/spacing)')
elseif strcmp(plotflag,'contourf')
contourf(tt,pseudofreq,abs(cwt_data),'edgecolor','none')
disp('Note: adapt number of contours if needed.')
end
%new_xtick_labels = cellfun(@(x) sprintf('%5.4f',x), num2cell(tt(get(gca,'XTick'))), 'uniformoutput', false);
%set(gca,'XTickLabel',new_xtick_labels)
%new_ytick_labels = cellfun(@(y) sprintf('%5.0f',y), num2cell(pseudofreq(get(gca,'YTick'))), 'uniformoutput', false);
%set(gca,'YTickLabel',new_ytick_labels,'YDir','reverse')
xlabel('Time (s)')
ylabel('Strouhal number')
xlim([tt(1) tt(end)])
ylim([pseudofreq(end) pseudofreq(1)])
%set(gca,'clim',[0 0.4*cwt_data_max])
end
%ondelette nous donne un diagramme de frequence en fonction de temps
%on a besoin le script 'monitorCWT_simple.m'
% % % % %n'a pas besoin de commander ''plot''
% % % % %faut calculer nbr de point, par example, dt=0.0084, St=0.23,
% % % % %period=1/St=4.34, nbr_point=period/dt=514, donc on prend 450:600
% % % % %'cmor1.5-1' c'est un caractere de ondelette, change pas ici!
[cwt_data pseudofreq]=monitorCWT_simple(monitor(:,3),burst,'cmor1.5-1',20:70,'contourf');
title('frequency modulation by wavelets');
colorbar;
%axis([0 80 1E-20 3]);
pseudofreq(end)
%---------------------------------------------------------------------------------------------------------------------------------------------------------
%----------------------------------------------AR-Autoregressive power spectral density estimate--------------------------------
%---------------------------------------------------------------------------------------------------------------------------------------------------------
clear all
clear y
close all
display('Please select the file')
monitor=uigetfile
monitor=importdata(monitor);
[burst]=burst_process(monitor(:,11)-mean(monitor(:,11)),2.2,monitor(:,3));
% %load('D:\bei\fft-AR\burst_process.m','tt','sig_burst','burst_HF','burst_BF')
% %load('/home1/emt2/dszubert/Simu/sa_2D_dense_monit/monitor_point/resultats/monitor_sig_defl_222.mat','tt','sig_defl')
sig=burst;
tt=monitor(:,3);
save_trace = 0;
font_size_label = 20;
font_name = 'Helvetica';
pAR = 200; % AR model order
% pARvect = 5:10:100;
% col = hsv(length(pARvect));
%sigtemp = figure('name','Signal temporel');
nfft = 1000; % FFT length calculation
fs = 1/(tt(2)-tt(1)); % Sampling frequency
dec_per = 4; % decalage periode
%per_ind = 12683;
per_ind =9172;
%deb = (12683*dec_per)+8000;
deb = 1;
sigtemp = figure(dec_per);
seg_length = 800;
nbr_seg = floor(per_ind/seg_length);
seg_length = floor(per_ind/nbr_seg); % longueur de segment ajustée pour prendre en compte le maximum de point
if seg_length < pAR
warning('seg_length < pAR')
end
col2 = hsv(nbr_seg);
%y(1,:) = sig_burst(deb:deb+per_ind,1)';
y(1,:) = sig(deb:deb+per_ind,1)';
ttloc = tt(deb:deb+per_ind)-tt(deb);
figure(1)
for k=1:nbr_seg
yseg(k,:) = y(1,(k-1)*seg_length+1:k*seg_length+1);
[Pxx(k,:),f(k,:)] = pyulear(yseg(k,:)-mean(yseg(k,:)),pAR,nfft,fs);
[PxxMax,indPxxMax] = max(Pxx(k,:));
fPxxMax(k) = f(k,indPxxMax);
loglog(f(k,:),Pxx(k,:))
pause
end
%% FFT
% Y = fft(y(1,:),nfft);
% Ymagn = abs(Y);
% fmax = fs/2;
% F1 = linspace(0,fmax,0.5*nfft);
%% Pwelch
%[Pxx2,F] = pwelch(sig_burst(:,1),hann(size(y,2)/2),65,nfft,fs);
%% Axes position
space = 0.01; % between graphs
margin = 0.1;
ax3 = [2*margin 2*margin 1-3*margin 0.75-2*margin-space];
ax1 = [2*margin 2*margin+ax3(4)+space ax3(3) 1-3*margin-space-ax3(4)];
%-ax2 = [margin margin 0.15 0.6];
%ax1 = [margin+ax2(3)+space margin+ax2(4)+space 1-2*margin-space-ax2(3) 1-margin-space-ax2(4)-0.08];
%ax3 = [margin+ax2(3)+space margin 1-2*margin-ax2(3)-space ax2(4)];
%% Figure finale
figure(sigtemp)
title('AR with filter Fc=2.2');
% Subplot signal
h1 = subplot(2,1,1);
plot(tt(deb:deb+per_ind)-tt(deb),y(1,:),'k','linewidth',1)
xlim([tt(deb) tt(deb+per_ind)]-tt(deb))
ymin = min(y(1,:));
ymax = max(y(1,:));
ampli = abs(ymin-ymax);
yliminf = ymin-0.1*abs(ampli);
ylimsup = ymax+0.1*abs(ampli);
ylim([yliminf ylimsup])
%title(['AR analysis - Windows length = ',num2str(seg_length), ' ; Signal length = ',num2str(size(y,2)),' ; AR order = ',num2str(pAR)],...
% 'FontSize',10,'FontName','FreeSans','FontWeight','bold')
%xlabel('Temps (s)')
ylabel('P_{norm}','FontWeight','bold','fontsize',font_size_label,'fontname',font_name)
hold on
% Mise en évidence des segments
for k=1:nbr_seg
plot([ttloc((k-1)*seg_length+1) ttloc((k-1)*seg_length+1)],[yliminf ylimsup],'k--') % Extremite inf
%plot([ttloc(k*seg_length+1) ttloc(k*seg_length+1)],[yliminf ylimsup],'--','color',col2(k,:)) % Extremite sup
end
hold off
set(h1,'position',ax1,'XTick',[],'fontsize',16,'fontname',font_name)
% h2 = subplot(2,2,3);
fPxxMaxamp = abs(max(fPxxMax)-min(fPxxMax));
fliminf = 1E-20; %min(fPxxMax)-0.4*abs(fPxxMaxamp);
if fliminf <= 0
fliminf = 1;
end
flimsup = 5; %max(fPxxMax)+0.4*abs(fPxxMaxamp);
% indtemp = find(F<fliminf);
% indliminf = indtemp(end);
% indtemp = find(F<flimsup);
% indlimsup = indtemp(end);
% %semilogx(Ymagn(indliminf:indlimsup),F(indliminf:indlimsup))
% semilogx(Pxx2(indliminf:indlimsup),F(indliminf:indlimsup))
%
% ylim([fliminf flimsup])
% set(gca,'XDir','reverse','Ygrid','on','YTickLabel',{});
% %xlabel('Spectre d''amplitude','FontWeight','bold')
% xlabel('PSD (Hanning)','FontWeight','bold')
% set(h2,'position',ax2)
h3 = subplot(2,1,2);
for k=1:nbr_seg
hold on
plot([ttloc((k-1)*seg_length+1) ttloc((k-1)*seg_length+1)],[fliminf flimsup],'k--') % Extremite inf
plot([ttloc((k-1)*seg_length+1) ttloc(k*seg_length+1)],[fPxxMax(k) fPxxMax(k)],'k','LineWidth',5.4)
hold off
end
ylim([fliminf flimsup])
xlim([tt(deb) tt(deb+per_ind)]-tt(deb))
xlabel('Time (s)','FontWeight','bold','fontsize',font_size_label,'fontname',font_name)
ylabel('Strouhal Number','FontWeight','bold','fontsize',font_size_label,'fontname',font_name)
set(gca,'YGrid','on');
set(h3,'position',ax3,'fontsize',16,'GridLineStyle','--','box','on','fontname',font_name)
%[~,title_handle] = suplabel(['AR analysis - Windows length = ',num2str(seg_length), ' ; Signal length = ',num2str(size(y,2)),' ; AR order = ',num2str(pAR)],'t');
%set(title_handle,'FontSize',10,'FontWeight','bold','fontname',font_name)
% A4 format
% xSize = 20.8; ySize = 29.5; % These are my size variables, width and height in cm
% xLeft = (21-xSize)/2; yTop = (30-ySize)/2; % Additional coordinates to center the figure on A4-paper
% set(gcf,'PaperUnits','centimeters') %This sets the units of the current figure (gcf = get current figure) on paper to centimeters
% set(gcf,'PaperPosition',[xLeft yTop xSize ySize])
if save_trace
print('-dpng','-r400',['/home1/emt2/dszubert/Signaux/Code_segmentation/AR/seg_length/',num2str(seg_length,'%06d'),'-',num2str(pAR,'%03d'),'.png'])
end
%clear yseg Pxx f fPxxMax
% figure(5)
% %
% % plot(Pxx2,F)
% % hold on
% % plot([min(Ymagn) max(Ymagn)],[F(indliminf) F(indliminf)],'--r')
% % plot([min(Ymagn) max(Ymagn)],[F(indlimsup) F(indlimsup)],'--r')
% % set(gca,'XScale','log','YScale','log','XDir','reverse','Ygrid','on');
% % hold off
% [cwt_data pseudofreq]=monitorCWT_simple(monitor(:,3),monitor(:,11)-mean(monitor(:,11)),'cmor1.5-1',20:200,'contourf');
% title('frequency modulation by wavelets');
% axis([0 80 pseudofreq(end) 3]);
% pseudofreq(end)
%---------------------------------------------------------------------------------------------------------------------------------------------------------
%-----------------------------------------------------------------------high-pass filter------------------------------------------------------------
%---------------------------------------------------------------------------------------------------------------------------------------------------------
function [burst] = burst_process(sig,Fc,tt)
close all
FlagPlot = 1;
FlagSavePlot = 0;
Fsize = 16; % taille titre graphs
% nbr_data = length(point);
% % frequence echantillonnage
% te = data_temp(2,1)-data_temp(1,1);
% fe = 1/te;
% % recherche frequence buffeting
% fenetre_freq = [70 90]; % freq buff = 78.8 Hz environ
% pmaxfft = monitormaxfft(k,point,data_fft,fenetre_freq);
% f_buff = sum(pmaxfft(:,3))/nbr_data;
% f_buff_std = std(pmaxfft(:,3));
% disp(['Frequence moyenne : ',num2str(f_buff), ' Hz (dispersion : ',num2str(f_buff_std),' Hz sur ',num2str(nbr_data),' valeurs)'])
% per_buff = 1/f_buff;
% clear pmaxfft
%% FFT
sig=sig-mean(sig);
yy=sig/sqrt(var(sig));
YY=fft(yy);
fe = 1/(tt(2)-tt(1));
dF=fe/length(yy);% length(yy)=length(a)
%F=0:dF:dF*length(yy)-dF;
%% EXTRACTION DES populations HF et BF BURSTS
disp('Extraction des populations HF et BF des bursts')
% load('monitor_sig.mat','tab_sig','tab_fft_sig')
% load('monitor_bursts_ifft.mat', 'tt', 'tab_sig_burst', 'F', 'tab_fft_burst', 'F_c_tremble', 'per_tremble')
% load monitor_sig_ifft tab_sig tab_fft_sig
% load monitor_bursts_ifft tt tab_sig_burst F tab_fft_burst F_c_tremble per_tremble
%load tremble_burst_ifft tab_sig tab_sig_burst tab_tremble tab_burst tab_fft_sig tab_fft_tremble tab_fft_burst t F rg_fin_tremble F_c_tremble
% DECOUPE (pour ne garder que la population HF BURSTS)
% Fc = 58/per_tremble; % 4.5726e+03
deb_burst = fix(Fc/dF)-3; % pour avoir var_i =0
fin_burst = length(sig)-deb_burst+3; % pour avoir var_i =0
disp([' separation bur./burst - frequence de coupure = ',num2str(Fc),' Hz'])
%rg_fin_HF = fin_burst; % 866;
% Ncc = size(sig,2);
% for t=1:Ncc;
YY_c = YY;
YY_c([1:deb_burst-1 fin_burst:length(YY_c)]) = 0;
x_c = ifft(YY_c);
burst = real(x_c);
fft_HF = YY_c;
% tab_burst_HF = [tab_burst_HF burst_HF];
% tab_fft_burst_HF = [tab_fft_burst_HF YY_c];
%
% Controle du caractere reel du signal
% var_i = var(imag(x_c))
% % Separation par DEFLATION pour ne garder que la population BF des BURSTS
%
% fft_burst = tab_fft_burst(:,t);
%
% YY_c = fft_burst-fft_HF;
% x_c = ifft(YY_c);
% q_c = real(x_c);
% burst_BF = q_c;
%
% tab_burst_BF = [tab_burst_BF burst_BF];
% tab_fft_burst_BF = [tab_fft_burst_BF YY_c];
% Controle du caractere reel du signal
% var_i = var(imag(x_c))
% end
% disp('Save monitor_bursts_ifft.mat...')
% save monitor_bursts_ifft tt tab_sig_burst tab_burst_HF tab_burst_BF F tab_fft_burst tab_fft_burst_HF tab_fft_burst_BF F_c_tremble Fc per_tremble
%save tremble_burst_ifft tab_sig tab_sig_burst tab_burst_HF tab_burst_BF tab_fft_sig tab_fft_tremble tab_fft_burst tab_fft_burst_HF tab_fft_burst_BF t F rg_fin_tremble F_c_tremble rg_fin_HF F_c_HF
%% TRACES SEPARES
disp('Traces separes')
if FlagPlot
Ncc = size(sig,2);
for t=1:Ncc;
% Le burst BF
% figure
% plot(tt,tab_burst_BF(:,t))
% v=axis;
% v(3)=1.2*v(3);
% v(4)=1.2*v(4);
% axis(v)
% title(['La partie BF du burst a F_c=', sprintf('%.0f Hz',F_c_HF),' au point ', num2str(t),' '],'fontsize',Fsize,'FontWeight','bold')
% xlabel('Time (s)','fontsize',Fsize*0.8,'FontWeight','bold')
% grid on
%
% if FlagSavePlot
% name = ['Sig_partie_BF_burst_m_ech_point_',num2str(t),'.jpg'];
% print('-djpeg',name)
% end
% Le burst HF
figure
plot(tt,burst(:,t))
title({'HF part of signal';['F_c = ', sprintf('%.2f ',Fc)]},'fontsize',Fsize,'FontWeight','bold')
xlabel('Time (s)','fontsize',Fsize*0.8,'FontWeight','bold')
ylabel('pressure variation (pa)')
grid on
if FlagSavePlot
name = ['Sig_partie_HF_burst_m_ech_point_',num2str(t),'.jpg'];
print('-djpeg',name)
end
end % for t
end % TRACES SEPARES