STEADY FLOW - TUTORIAL 2 FOR MSc 1st year in Paul Sabatier University

This script has been designed to study the flow around a cylinder
at low Reynolds number [1-40].
We study how this flow is characterised by analysing main properties
of flow quantities such as flow field, pressure, vorticity and
shear stress.
The script performs the following calculations :
# Generation of a mesh
# Computation of the stationary solution at Re = [1-40]
# Display flow fields

Contents

First : initialization of StabFem tools

close all;
addpath('../../SOURCES_MATLAB');
SF_RecoverDataBase('WORK_CYLINDER');
SF_Start('verbosity', 3,'ffdatadir', './WORK_CYLINDER/');
Directory WORK_CYLINDER is aready present
NOTICE  - Initialization already done
NOTICE  - Verbosity set to 3 : Notice messages
NOTICE  - Database directory :./WORK_CYLINDER/
NOTICE  - This directory already exists
NOTICE  - Please type "SF_Status" to see what is available in this dir
NOTICE  -     or type "SF_core_arborescence('cleanall')" to erase any previous data

Check what is available in the dataset

SF_Status

% Here we see that Re = 40 is the field number 7
% Hence we load it this way :

 bf = SF_Load('BASEFLOWS',7);
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    ./WORK_CYLINDER/
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK_CYLINDER/MESHES :
     (list of meshes previously created/adapted ; couples of .msh/.ff2m files )
 
Index | Name              | generation mode | Date                 | Nv      
1     | FFMESH_000001.msh | initial         | 23-mars-2023 11:44:59 | 2192    
2     | FFMESH_000002.msh | adapted         | 23-mars-2023 11:45:44 | 1986    
3     | FFMESH_000003.msh | adapted         | 23-mars-2023 11:45:52 | 3340    
4     | FFMESH_000004.msh | mirror          | 23-mars-2023 11:46:00 | 6544    
#################################################################################
 
     (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         | Fx        
 1     | FFDATA_000003.txt | BaseFlow     | 23-mars-2023 11:45:46 | FFMESH_000002.msh | 60         | 0.65237   
 2     | FFDATA_000004.txt | BaseFlow     | 23-mars-2023 11:45:55 | FFMESH_000003.msh | 60         | 0.65237   
 3     | FFDATA_000011.txt | BaseFlow     | 23-mars-2023 11:46:07 | FFMESH_000004.msh | 60         | 0.63298   
 
#################################################################################
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK_CYLINDER/BASEFLOWS
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Fx        
 1     | FFDATA_000001.txt | BaseFlow     | 23-mars-2023 11:45:14 | FFMESH_000001.msh | 1          | 5.5507    
 2     | FFDATA_000002.txt | BaseFlow     | 23-mars-2023 11:45:25 | FFMESH_000001.msh | 10         | 1.4307    
 3     | FFDATA_000003.txt | BaseFlow     | 23-mars-2023 11:45:38 | FFMESH_000001.msh | 60         | 0.65237   
 4     | FFDATA_000004.txt | BaseFlow     | 23-mars-2023 11:46:44 | FFMESH_000004.msh | 60         | 0.64404   
 5     | FFDATA_000005.txt | BaseFlow     | 23-mars-2023 11:47:29 | FFMESH_000004.msh | 40         | 0.75381   
 6     | FFDATA_000006.txt | BaseFlow     | 23-mars-2023 11:48:08 | FFMESH_000004.msh | 10         | 1.3906    
 7     | FFDATA_000007.txt | BaseFlow     | 23-mars-2023 11:48:51 | FFMESH_000004.msh | 1          | 5.3375    
 8     | FFDATA_000008.txt | BaseFlow     | 23-mars-2023 11:49:38 | FFMESH_000004.msh | 40         | 0.75381   
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK_CYLINDER/MESHBF
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Fx        
 1     | FFDATA_000001.txt | BaseFlow     | 23-mars-2023 11:45:13 | FFMESH_000001.msh | 1          | 5.5507    
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK_CYLINDER/MISC
     (couples of .txt/.ff2m files )
   NB this directory may contain various secondary files produced by StabFem,
      such as flowfields projected after adaptation but not recomputed, adaptation masks, etc... 
 
 Index | Name              | Type         | Date                 | Mesh file        
 1     | FFDATA_000001.txt | MASK         | 23-mars-2023 11:45:42 | FFMESH_000001.msh
 2     | FFDATA_000002.txt | MASK         | 23-mars-2023 11:45:50 | FFMESH_000002.msh
 
 
#################################################################################

ans = 

  struct with fields:

       MESHES: [1x4 struct]
    BASEFLOWS: [1x8 struct]
       MESHBF: [1x1 struct]
         MISC: [1x2 struct]

NOTICE  - Loading baseflow number 7

plot Ux, Uy, omega, P

Il montre le dernier baseflow. Ici vous pouvez afficher vos champs.

figure();
subplot(2,2,1); SF_Plot(bf,'ux','xlim',[-3 5],'ylim',[-3 3],'title','u_x'); hold on;
SF_Plot(bf,'psi','Contour','only','CLevels',[-0.05:0.01:0.05],'CColor','w','xlim',[-3 5],'ylim',[-3 3]);
subplot(2,2,2); SF_Plot(bf,'uy','xlim',[-3 5],'ylim',[-3 3],'title','u_y'); hold on;
SF_Plot(bf,'psi','Contour','only','CLevels',[-0.05:0.01:0.05],'CColor','k','xlim',[-3 5],'ylim',[-3 3]);
subplot(2,2,3); SF_Plot(bf,'vort','Contour','on','CLevels',[-1,-0.5,-0.02,0.02,0.5,1],'xlim',[-3 5],'ylim',[-3 3],'title','\omega','colorrange',[-1,1]);
subplot(2,2,4); SF_Plot(bf,'p','Contour','on','xlim',[-3 5],'ylim',[-3 3],'title','p');
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Quiver plot

bf.Normu = sqrt(bf.ux.^2+bf.uy.^2);
% figure ; SF_Plot(bf,'Normu')

figure;
SF_Plot(bf,'Normu','xlim',[-1 5],'ylim',[-3 3],'title','u_x');
hold on;
SF_Plot(bf,'psi','Contour','only','CLevels',[-0.05:0.01:0.05],'CColor','y','xlim',[-3 5],'ylim',[-3 3]);
SF_Plot(bf,{'ux','uy'},'xlim',[-1 5],'ylim',[-3 3],'FColor','w');
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

If you want to plot the norm of the velocity

% Velocity plots

% tout le long du cylindre (RMin = 0.5). Vous pouvez modifier la valeur de
% RMin pour ainsi analyser d'un point de vue local l'écoulement autour le
% cylindre
RMax = 50; RMin = 0.5; % RMin et RMax sont modifiables. Rappelez vous RMax <= min(xmax,ymax)
R = [RMin:.01:RMax];
thetaList = [pi/12:pi/12:pi/2];
nx = size(thetaList); nx = nx(2);
ny = size(R); ny = ny(2);

UBLVec = zeros(nx,ny);
VBLVec = zeros(nx,ny);
UrBLVec = zeros(nx,ny);
UthBLVec = zeros(nx,ny);

for index=[1:nx]
    theta = thetaList(index);
    X=R.*cos(theta); Y=R.*sin(theta);
    V=SF_ExtractData(bf,'uy',X,Y); VBLVec(index,:) = V;
    U=SF_ExtractData(bf,'ux',X,Y); UBLVec(index,:) = U;
    UrBL = U.*cos(theta) + V.*sin(theta); UrBLVec(index,:) = UrBL;
    UthBL = -U.*sin(theta) + V.*cos(theta); UthBLVec(index,:) = UthBL;
end

% U_theta - Couche limite du cylindre (coord polaires)
figure;
hold on;
for index=[1:nx]
    V = UthBLVec(index,:);
    legendName = ['$\theta = ',num2str(index),'\pi/12 $'];
    plot(V,R,'o','DisplayName',legendName);
    xlabel('$U_{\theta}(x,y)$','Interpreter','latex');
    ylabel('$r$','Interpreter','latex');
    title('$U_{\theta}(r,\theta)$', 'Interpreter','latex');
end
leg = legend;
set(leg,'Interpreter','latex','fontsize',24);
ylim([0,4]);
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Calcul de la contrainte parietale

theta = linspace(0,2*pi,200);
Ycircle = 0.5001*sin(theta); Xcircle = 0.5001*cos(theta);

tauxxwall = SF_ExtractData(bf,'tauxx',Xcircle,Ycircle);
tauxywall = SF_ExtractData(bf,'tauxy',Xcircle,Ycircle);
tauyywall = SF_ExtractData(bf,'tauyy',Xcircle,Ycircle);
S = sin(theta); C = cos(theta);
tauWallN =  C.*(tauxxwall.*C+tauxywall.*S) + S.*(tauxywall.*C+tauyywall.*S);
tauWallT = - S.*(tauxxwall.*C+tauxywall.*S) + C.*(tauxywall.*C+tauyywall.*S);
pWall = SF_ExtractData(bf,'p',Xcircle,Ycircle);

figure;
plot(theta,tauWallT,'b^');
hold on;
plot(theta,tauWallN,'r^');
xlabel('\theta'); ylabel('\tau_{w}'); title('shear stress $\tau(r=R,\theta)$ along the surface', 'Interpreter','latex');
legend({'$\tau_t$','$\tau_n$'},'Interpreter','latex');


% Calcul de la traînée
Cx = 2*trapz(theta, 0.5*((tauWallN.*C - tauWallT.*S)/bf.Re - pWall.*C))

%%[[PUBLISH]]


%
fclose(fopen('SCRIPT_Re40.success','w+'));
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Cx =

   10.7576