LINEAR Stability ANALYSIS of the wake of a SPRING-MOUNTED CYLINDER

this scripts demonstrates how to use StabFem to study the instabilities
in a Vortex-Induced Vibrations (VIV) problem.
The script reproduces figures 2c, 4a and 5a from the paper by D. Sabino et al.

Contents

First we set a few global variables needed by StabFem

addpath([fileparts(fileparts(pwd)),'/SOURCES_MATLAB']);
SF_Start('verbosity',4);
close all;
NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Compute or import the base flow

bf = SF_Load('lastadapted');

if isempty(bf)
    bf = SmartMesh_Cylinder;
end
 
mesh adaptation  : type S
 Adapted mesh has been generated ; number of vertices = 2218

Eigenvalue computation : plot one mode

bf = SF_BaseFlow(bf,'Re',60);
mstar = 20; Ustar = 3;shiftFM = 0.04689 + 0.74874i;
M = mstar*pi/4; K = pi^3*mstar/Ustar^2;
[ev,em] = SF_Stability(bf,'shift',shiftFM,'nev',1,'type','D','STIFFNESS',K,'MASS',M,'DAMPING',0);
figure ; SF_Plot(em(1),'ux','xlim',[-2 4],'ylim', [0 2],'colorrange','cropcentered');

Loop on eigenvalue computation for Mstar = 20 (figure 4a)

starting points

shiftFM = 0.04689 + 0.74874i;
shiftEM = -0.0188 + 2.0100i;

% computing branch "FM" for a range of Ustar
Ustar_tab = [3:.2:11];
K = pi^3*mstar/Ustar_tab(1)^2;
ev = SF_Stability(bf,'shift',shiftFM,'nev',1,'type','D','STIFFNESS',K,'MASS',M,'DAMPING',0);

for k = 1:length(Ustar_tab)
  Ustar = Ustar_tab(k);
  K = pi^3*mstar/Ustar^2;
  ev = SF_Stability(bf,'shift','cont','nev',1,'type','D','STIFFNESS',K,'MASS',M,'DAMPING',0);
  FM_tab(k) = ev;
end

% computing branch "EM" for a range of Ustar
Ustar_tab = [3:.2:11];
K = pi^3*mstar/Ustar_tab(1)^2;
ev = SF_Stability(bf,'shift',shiftEM,'nev',1,'type','D','STIFFNESS',K,'MASS',M,'DAMPING',0);

for k = 1:length(Ustar_tab)
  Ustar = Ustar_tab(k);
  K = pi^3*mstar/Ustar^2;
  ev = SF_Stability(bf,'shift','cont','nev',1,'type','D','STIFFNESS',K,'MASS',M,'DAMPING',0);
  EM_tab(k) = ev;
end

% plotting both these branches
figure(30);
subplot(2,1,1);
plot(Ustar_tab,real(FM_tab),'b');hold on;
plot(Ustar_tab,real(EM_tab),'r');
xlabel('U*'); ylabel('\lambda_r');
legend('FM','EM');
subplot(2,1,2);
plot(Ustar_tab,imag(FM_tab)/(2*pi),'b');hold on;
plot(Ustar_tab,imag(EM_tab)/(2*pi),'r');
xlabel('U*'); ylabel('St');
legend('FM','EM');

Impedance computation and plot for Re = 35

bf = SF_BaseFlow(bf,'Re',35);
fo = SF_LinearForced(bf,'omega',[.01:.01:1],'plot','yes');

Impedance computation and plot for Re = 60

bf = SF_BaseFlow(bf,'Re',60);
fo = SF_LinearForced(bf,'omega',[.02:.02:2],'plot','yes');

% Comparison between eigenvalue and impedance-based approximation for Re =
% 60, m* = 20 (figure 5a)
mstar = 20; gamma = 0;
fo = Uasympt(fo,mstar,gamma); % See this function
figure(30);
subplot(2,1,1);hold on;plot(fo.Ustar,fo.omegai,'k--');xlim([3 11]);legend('off');
subplot(2,1,2);hold on;plot(fo.Ustar,fo.omega/(2*pi),'k--');xlim([3 11]);
legend('Fluid mode','Elastic mode','Impedance-based prediction');
Warning: Imaginary parts of complex X and/or Y arguments ignored. 

[[PUBLISH]] (this tag is to enable automatic publication as html ; don't touch this line)

%
fclose(fopen('VIV.success','w+'));