Parametric study of the flow around a confined cylinder (part 1)

This example script shows how to perform a parametric study of a problem. We consider a cylinder of diameter 1 placed in a plane chanel of width H. We want to track the thresholds as function of Re and H.

This program performs the parametric study and stores the data in a database folder.

The postprocessing is done in a second program : see ParametricStudy_PostProcess.m

Note the tag [[WORK]] which is needed to tell the server to do a backup of the database folder 'WORK_ConfinedCyl'. This backup will be reopened by the posprocessing program.

Contents

Chapter 0 : Initialization

addpath('../../SOURCES_MATLAB/');
SF_Start('verbosity',2);
SF_DataBase('create','WORK_ConfinedCyl');

NOTICE  - Initialization already done

Outer Loop over H

shiftStart = 0.14471+1.2152i;

for H = [10 7.5 5 4 3 2.5 2]

  H

H =

    10

H =

    7.5000

H =

     5

H =

     4

H =

     3

H =

    2.5000

H =

     2

Create mesh, compute baseflow and adapt to one eigenmode (adjoint) for Re = 60

  ffmesh =  SF_Mesh('Mesh_Cylinder.edp','Options',{'Xmin',-10,'Xmax',30,'Ymax',H/2,'labelsides',11});
  bf = SF_BaseFlow(ffmesh,'Solver','Newton_2D.edp','Options',{'Re',1});
  bf = SF_Adapt(bf);
  bf = SF_BaseFlow(bf,'Options',{'Re',10});
  bf = SF_BaseFlow(bf,'Options',{'Re',60});
  [ev,em]= SF_Stability(bf,'shift',shiftStart,'nev',5,'type','A','solver','Stab_2D.edp','Stats',false);
  % NB here option 'Stats',false is to avoid writing this eigenvalue in the
  % statistics, hence it won't appear on the final curves
  bf = SF_Adapt(bf,em(1));
  bf = SF_BaseFlow(bf,'Options',{'Re',60});
  [ev,em]= SF_Stability(bf,'shift',conj(ev(1)),'nev',1,'type','D','solver','Stab_2D.edp','Stats',false);
  shiftStart = ev(1); % shift for next eigenmode needed for adapting mesh

WARNING -  Folder WORK_ConfinedCyl/MESHBF not found : creating it !

plotting the baseflow and the eigenmode (for H = 10 )

  if (H == 10 )
    figure; SF_Plot(bf,'ux');hold on; SF_Plot(bf,'psi','contour','only','xlim',[-2 5],'ylim',[ 0 5])
    figure; SF_Plot(em,'ux','xlim',[-2 5],'ylim',[ 0 5],'colorrange','cropcenter')
  end

%plotoptions = {'xlim' [-2 5]};
%[ev,em]= SF_Stability(bf,'shift','cont','nev',10,'Stats',false,'Plotspectrum',true,'PlotModeOptions',plotoptions);

ans = 

  2x1 Contour array:

  Contour
  Contour


ans = 

  2x1 graphics array:

  ColorBar
  Patch

Inner Loop over Re

  tt = 'off';shift = shiftStart; % disable threshold tracking for first computation
  Recguess = 25+250*H^(-1.5); % guess coming from interpolation of a few preliminary computations
  for Re = Recguess*[1.5:-.125:.75]
     bf = SF_BaseFlow(bf,'Options',{'Re',Re});
     [ev,em]= SF_Stability(bf,'shift','cont','nev',10,'Threshold',tt);
     tt = 'single';shift='cont'; % enable threshold tracking
  end

WARNING -  Folder WORK_ConfinedCyl/NEARLYMARGINALMODES not found : creating it !

end

Plot results

sfs = SF_Status('QUIET')

sfs = 

  struct with fields:

                 MESHES: [1x21 struct]
              BASEFLOWS: [1x77 struct]
             EIGENMODES: [1x532 struct]
            EIGENVALUES: [1x1 struct]
             THRESHOLDS: [1x1 struct]
                 MESHBF: [1x7 struct]
                   MISC: [1x14 struct]
    NEARLYMARGINALMODES: [1x6 struct]

First all eigenvalues detected as function of Re

figure ;
subplot(2,1,1); plot(sfs.EIGENVALUES.Re,real(sfs.EIGENVALUES.lambda),'x',sfs.THRESHOLDS.Re,real(sfs.THRESHOLDS.lambda),'ro')
xlabel('Re');ylabel('\lambda_r');ylim([-.1 .4])
subplot(2,1,2); plot(sfs.EIGENVALUES.Re,imag(sfs.EIGENVALUES.lambda).*(real(sfs.EIGENVALUES.lambda)>-0.1),'x',sfs.THRESHOLDS.Re,imag(sfs.THRESHOLDS.lambda),'ro')
xlabel('Re');ylabel('\lambda_i');ylim([.75 3.5])

Plot Re_c as function of H/D

figure;
plot(sfs.THRESHOLDS.Re,sfs.THRESHOLDS.H);xlabel('Re_c'); ylabel('H/D');

Note on geometry and implementation of boundary conditions

We use half domain so dimensions is y in [0-H/2] ; x in [-10,30]

The inlet profile is plane Poiseuille with mean velocity Um = 1. Hence the inlet velocity profile is Uinlet = 3/2*(1-(2*y/H)^2)

This definition is done through macro Uinlet defined in SF_Custom.idp. SF_Custom.idp also contains redefinitions of some macros : SFWriteMesh : definition of H as a metadata for the mesh SFWriteBaseFlow : note that we have to modify boundary conditions for streamfunction psi to fit with inlet condition.

% [[PUBLISH]]
% [[WORK]]

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

Published with MATLAB® R2022b