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