Computation of the flow around a rotating cylinder using pseudo-arclength continuation

This program reproduces figure 3(c) from Sierra et al (JFM, 2020)

Contents

Initialise

close all;
addpath('../../SOURCES_MATLAB');
SF_Start('verbosity', 4);
SF_DataBase('create','./WORK');
SF_core_setopt('ErrorIfDiverge',false); % option recommended if netwon divergence is correctly caught.
NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Generation of initial solution and mesh

ffmesh = SF_Mesh('Mesh_Cylinder_FullDomain.edp','Params',[-40 80 40],'problemtype','2d');
bf=SF_BaseFlow(ffmesh,'Re',1,'Omegax',4.5);
for Re = [10, 60, 100 170]
    bf=SF_BaseFlow(bf,'Re',Re,'Omegax',4.5);
end
bf = SF_Adapt(bf,'anisomax',10,'Hmax',10);
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/MESHBF not found : creating it !

Arclength continuation

alpha = bf.Omegax;
alphaMax = 5.5; stepMax = 0.1;
step = stepMax;
FxList = []; FyList = []; AlphaList = [];

while ( bf.Omegax < alphaMax)
    bf=SF_BFContinue(bf,'step',step);
    % NB in this example we use the solver 'ArcLengthContinuation2D.edp'
    % which is the default one for this class of problems.
    while(bf.iter < 0)
        SF_core_log('w','Not converged. Arc-length step reduced by a factor of 2');
        step = 0.5*step
        bf=SF_BFContinue(bf,'step',step);
        if (abs(step)<stepMax/100)
            warning('step too small')
            continue;
        end
    end
    step = sign(step)*min(1.4*abs(step),stepMax);
    FxList = [FxList, bf.Fx];
    FyList = [FyList, bf.Fy];
    AlphaList = [AlphaList, bf.Omegax];
end
ERROR   - SF_core_freefem: Error while using FreeFem++  ; file = ArcLengthContinuation2D.edp ; error code = 1 : Syntax error in your .edp program
WARNING - An error is encountered while publishing script. Launching SF_Status for diagnostics
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    ./WORK/
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MESHES :
     (list of meshes previously created/adapted ; couples of .msh/.ff2m files )
 
Index | Name              | generation mode | Date                 | Nv      
1     | FFMESH_000001.msh | initial         | 08-sept.-2021 01:05:35 | 3908    
2     | FFMESH_000002.msh | adapted         | 08-sept.-2021 01:06:35 | 2173    
#################################################################################
 
     (list of base flows associated to newly computed meshes ; couples of .txt/.ff2m files )
     NB : these base flows are most likely simply projected and not recomputed, it is not recommended to use them
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Omegax     | Fx         | Fy         | beta      
 1     | FFDATA_000003.txt | BaseFlow     | 08-sept.-2021 01:06:36 | FFMESH_000002.msh | 170        | 4.5        | 0.056888   | 11.0601    | 0.28444   
 
#################################################################################
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/BASEFLOWS
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Omegax     | Fx         | Fy         | beta      
 1     | FFDATA_000001.txt | BaseFlow     | 08-sept.-2021 01:05:48 | FFMESH_000001.msh | 1          | 4.5        | 4.7358     | 7.6589     | 23.6789   
 2     | FFDATA_000002.txt | BaseFlow     | 08-sept.-2021 01:05:59 | FFMESH_000001.msh | 10         | 4.5        | 0.18368    | 8.7724     | 0.91841   
 3     | FFDATA_000003.txt | BaseFlow     | 08-sept.-2021 01:06:12 | FFMESH_000001.msh | 60         | 4.5        | -0.020816  | 10.2893    | -0.10408  
 4     | FFDATA_000004.txt | BaseFlow     | 08-sept.-2021 01:06:22 | FFMESH_000001.msh | 100        | 4.5        | 0.016006   | 10.6732    | 0.080031  
 5     | FFDATA_000005.txt | BaseFlow     | 08-sept.-2021 01:06:33 | FFMESH_000001.msh | 170        | 4.5        | 0.056888   | 11.0601    | 0.28444   
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MESHBF
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Omegax     | Fx         | Fy         | beta      
 1     | FFDATA_000001.txt | BaseFlow     | 08-sept.-2021 01:05:48 | FFMESH_000001.msh | 1          | 4.5        | 4.7358     | 7.6589     | 23.6789   
 
#################################################################################
Error using SF_core_log (line 124)
SF_core_freefem: Error while using FreeFem++  ; file = ArcLengthContinuation2D.edp ; error code = 1 : Syntax error in your .edp program

Error in SF_core_freefem (line 409)
     SF_core_log(errortype, ['SF_core_freefem: Error while using FreeFem++  ; file = ', cmd,' ; error code = ',num2str(s) ,' : ' ,errormsg]);

Error in SF_BFContinue (line 154)
value = SF_core_freefem(ffsolver,'parameters',ffparams,'arguments',ffarguments);

Error in ArcLength (line 32)
    bf=SF_BFContinue(bf,'step',step);

Plot

Fy

h = figure;
subplot(1,2,1);
plot(AlphaList,FxList,'ok-','linewidth',2)
set(gca,'TickLabelInterpreter','latex','fontsize',24);
xlabel('$\alpha$','interpreter','Latex');
ylabel('$F_x$','interpreter','Latex');

set(gcf, 'Color', 'w');

subplot(1,2,2);
plot(AlphaList,FyList,'ok-','linewidth',2)
set(gca,'TickLabelInterpreter','latex','fontsize',24);
xlabel('$\alpha$','interpreter','Latex');
ylabel('$F_y$','interpreter','Latex');
set(gcf, 'Color', 'w');


% Resize automatically for saving
fig = gcf;
fig.PaperPositionMode = 'auto'
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
set(fig,'Units','Inches');
print(fig,'FyAlpha','-dpdf','-bestfit');

Visualization

index = 20;
bf = SF_Load('BASEFLOWS',index);
bf.Omegax
figure;
Xmin = -2; Xmax = 2; Ymin = -2; Ymax = 4;
SF_Plot(bf,'vort','xlim',[Xmin,Xmax],'ylim',[Ymin,Ymax]);
hold on;
SF_Plot(bf,'psi','contour','only','CLevels',[-10.0:0.05:10],...
    'xlim',[Xmin,Xmax],'ylim',[Ymin,Ymax]);

Stability calculation

[ev,em] = SF_Stability(bf,'shift',0.1,'nev',5,'plotspectrum',true);

Summary

SF_Status

% [[PUBLISH]]

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