Computation of base flow and DNS simulation for the 2D incompressible flow around a cylinder.

This program is used to recompute all data provided in the practical work at UPS (M1 Mecanique)

In this tutorial script: 1/ we bluild a well-suited mesh 2/ we compute base-flows (steady solutions) for sample values of Re 3/ we perform a DNS of the wake of a cylinder for Re = 60, allowing to observe the development and saturation of the well-known Von Karman vortex street.

The FreeFem++ code TimeStepper_2D.edp is adapted from the one used in Jallas, Marquet & Fabre (PRE, 2017).

Chapter 0 : initialization
Chapter 1 : mesh
Chapter 2 : steady flow
Plot the solution for Re = 40
Intermezzo : compute a few potential solutions
Chapter 3
Then we Launch a DNS
second run over one period to generate snapshots ;
Chapter 4 : plotting the results and generating a movie
Summary

Chapter 0 : initialization

clear all ; close all;
addpath('../../SOURCES_MATLAB');
SF_Start('verbosity', 4);
SF_DataBase('create', './WORK_CYLINDER/');


% NB change verbosity to 4 to follow the simulation while it's running
%  or 2 to supress output from freefem++

 Warning  : SF_core_start was not previously lanched ; launching right now
NOTICE  - Verbosity set to 3 : Notice messages
Startup 01: Checking MATLABPATH
Startup 02: Checking Octave Path
Startup 03: Loading live options
NOTICE  - Verbosity set to 3 : Notice messages
NOTICE  - Git is available on this platform.
NOTICE  - Successfully detected FreeFem++-mpi
Startup 04: Loading static options
NOTICE  - SF_core_opts/read: stabfem.opts could not be found
Startup 05: Post-startup tasks
NOTICE  - Detected FreeFem++-mpi version 4.9
Startup 06: FreeFem tests 
NOTICE  - librairy SuperLu is available
NOTICE  - librairy MUMPS is available
NOTICE  - librairy PETSc is available
NOTICE  - librairy SLEPc is available
NOTICE  - librairy SLEPc-complex is available
NOTICE  - Using SLEPc-complex as eigensolver
WARNING - librairy ARPACK NOT available
NOTICE  - Elementary Freefem test 1 (import data) passed
NOTICE  - Elementary Freefem test 2 (piping parameters) passed
NOTICE  -  Startup ENDED ; StabFem is now ready for use !
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Chapter 1 : mesh

We use the same procedure as in CYLINDER_LINEAR.m. First we build an initial mesh and progressively increase the Reynolds :

ffmesh = SF_Mesh('Mesh_Cylinder.edp','Params',[-40 80 40],'problemtype','2D','symmetry','S');
bf=SF_BaseFlow(ffmesh,'Re',1);
bf=SF_BaseFlow(bf,'Re',10);
bf=SF_BaseFlow(bf,'Re',60);
Mask = SF_Mask(bf.mesh,[-2 10 0 2 .15]);
bf = SF_Adapt(bf,Mask,'Hmax',5);
Mask = SF_Mask(bf.mesh,[-2 10 0 2 .15]);
bf = SF_Adapt(bf,Mask,'Hmax',5);
bf = SF_Mirror(bf);

WARNING -  Your program uses depreciated syntax 'Params' to pass parameters 
WARNING -  It is advised to switch to new method using 'Options', and modify your Freefem solvers to detect parameters using 'getARGV'  
WARNING -  Folder ./WORK_CYLINDER/MESHBF not found : creating it !

Chapter 2 : steady flow

Now computes the base flow for a few values of Re needed for the TP

bf = SF_BaseFlow(bf,'Re',60);
bf = SF_BaseFlow(bf,'Re',40);
bf = SF_BaseFlow(bf,'Re',10);
bf = SF_BaseFlow(bf,'Re',1);

Plot the solution for Re = 40

bf = SF_BaseFlow(bf,'Re',40);

SF_Plot(bf,'ux','xlim',[-2 10],'ylim',[-3 3 ],'colorrange','cropminmax','boundary','on','bdlabels',2,'title','Steady flow for Re = 40');
hold on; SF_Plot(bf,'psi','contour','only','clevels',[-2 -1.5 -1 -.5  -.2 -.1 -0.05 -.02 0 0.02 0.05 0.1 .2 .5 1 1.5 2],'xlim',[-2 10],'ylim',[-3 3]);

Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Intermezzo : compute a few potential solutions

mesh = SF_Load('MESH','last');
FlowP0 = SF_Launch('PotentialFlow.edp','Options',{'U',1,'V',0,'Gamma',0},'mesh',mesh,...
                  'DataFile','PotentialFlow.txt','Store','POTENTIALFLOWS','LoadFiles',false);
FlowP1 = SF_Launch('PotentialFlow.edp','Options',{'U',1,'V',0,'Gamma',pi},'mesh',mesh,...
                  'DataFile','PotentialFlow.txt','Store','POTENTIALFLOWS');
FlowP2 = SF_Launch('PotentialFlow.edp','Options',{'U',1,'V',0,'Gamma',2*pi},'mesh',mesh,...
                  'DataFile','PotentialFlow.txt','Store','POTENTIALFLOWS');
FlowP3 = SF_Launch('PotentialFlow.edp','Options',{'U',1,'V',0,'Gamma',3*pi},'mesh',mesh,...
                  'DataFile','PotentialFlow.txt','Store','POTENTIALFLOWS');

figure; SF_Plot(FlowP3,'p'); hold on; SF_Plot(FlowP3,{'ux','uy'},'xlim',[-2 2],'ylim',[-2 2])

WARNING - Obsolete usage of SF_Load : must specify a folder

Error using SF_Launch
'LoadFiles' is not a recognized parameter. For a list of valid name-value pair
arguments, see the documentation for this function.

Error in SCRIPT_GENERATE_DATA (line 60)
FlowP0 = SF_Launch('PotentialFlow.edp','Options',{'U',1,'V',0,'Gamma',0},'mesh',mesh,...

Chapter 3

We first generate an initial condition for the DNS at Re = 60 using linear stability results

bf = SF_BaseFlow(bf,'Re',60);
[ev,em] = SF_Stability(bf,'shift',0.04+0.76i,'sym','N','nev',1); % compute the eigenmode.
startfield = SF_Add({bf,em(1)},'coefs',[1 0.01],'NameParams',['Re'],'Params',[bf.Re]); % creates startfield = bf+0.01*em

Then we Launch a DNS

We use the driver SF_DNS.m which launches the Freefem++ solver TimeStepper_2D.edp

% first run over a long period without snapshots

Niter = 8000;   % total number of time steps
iout = 8000;   % will generate snapshots each iout time steps
dt = 0.02;   % time step
DNSOptions = {'Re',60,'Niter',Niter,'dt',dt,'iout',iout};

%[DNSstatslong,DNSfieldslong] =SF_DNS(startfield,'Re',60,'Nit',Nit,'dt',dt,'iout',iout); %first run on 8000 time steps

folder = SF_TS_Launch('TimeStepper_2D.edp','Init',startfield,'Options',DNSOptions,'wait',true);
[DNSstatslong,DNSfieldslong] = SF_TS_Check(folder); % imports results

second run over one period to generate snapshots ;

dt = 7.4/400; % to simulate one period in 400 time steps
iout = 40;
Niter=400;
DNSOptions = {'Re',60,'Niter',Niter,'dt',dt,'iout',iout};
startfield = DNSfieldslong(end);
%[DNSstats,DNSfields] = SF_DNS(DNSfieldslong(end),'Re',60,'Nit',Nit,'dt',dt,'iout',40); %second run on approximately one period
folder = SF_TS_Launch('TimeStepper_2D.edp','Init',startfield,'Options',DNSOptions,'wait',true);
[DNSstats,DNSfields] = SF_TS_Check(folder); % imports results

Chapter 4 : plotting the results and generating a movie

% Here is how to generate a movie
h = figure;
mkdir('html');
filename = 'html/DNS_Cylinder_Re60.gif';
for i=1:length(DNSfields)
    SF_Plot(DNSfields(i),'vort','xlim',[-2 10],'ylim',[-3 3 ],'colorrange',[-2 2],...
        'title',['t = ',num2str(DNSfields(i).t)],'boundary','on','bdlabels',2);
    hold on; SF_Plot(DNSfields(i),'psi','contour','only','clevels',[-2 -1.5 -1 -.5  -.2 -.1 -0.05 -.02 0 0.02 0.05 0.1 .2 .5 1 1.5 2],'xlim',[-2 10],'ylim',[-3 3]);
    hold off;
    pause(0.01);
    frame = getframe(h);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if (i==1)
       imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.1);
       set(gca,'nextplot','replacechildren');
    else
       imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.1);
    end
end

Here is the movie

We now plot the lift force as function of time

figure(15);
subplot(2,1,1);
plot(DNSstatslong.t,DNSstatslong.Fy);title('Lift force as function of time');
xlabel('t');ylabel('Fy');
subplot(2,1,2);
plot(DNSstatslong.t,DNSstatslong.Fx);title('Drag force as function of time');
xlabel('t');ylabel('Fx');

Now we plot the pressure and vorticity along the surface

alpha = linspace(-pi,pi,201);
Xsurf = .501*cos(alpha);
Ysurf = .501*sin(alpha);
Psurf = SF_ExtractData(DNSfields(end),'p',Xsurf,Ysurf);
Omegasurf = SF_ExtractData(DNSfields(end),'vort',Xsurf,Ysurf);
figure(16);
subplot(2,1,1);title('Pressure along the surface P(r=a,\theta) at final time step')
plot(alpha,Psurf);
xlabel('\theta');ylabel('P(r=a,\theta)');
subplot(2,1,2);title('Vorticity along the surface \omega(r=a,\theta) at final time step')
plot(alpha,Omegasurf);
xlabel('\theta');ylabel('\omega_z(r=a,\theta)');

Summary

SF_Status


% [[PUBLISH]]
% [[WORK]]
%
fclose(fopen('SCRIPT_GENERATE_DATA.success','w+'));