Computation of the potential flow around a 2D airfoil with FreeFem/StabFem

This example is a demonstration of how to compute a potential flow or a viscous flow around a wing profile (NACA1408).

The source code of this program is Example_NACA1408.m

First initialize StabFem:
Generation of a mesh
Compute Potential flow for alpha = 0 ° ; Gamma = 0 (Kutta condition not verified)
Computing flow for alpha = 0 degrees with right circulation (Kutta condition verified)
VISCOUS FLOWS
Compare the Lift and drag forces in potential and viscous cases
SUMMARY

First initialize StabFem:

close all;
addpath('../../SOURCES_MATLAB');
SF_Start('verbosity',4,'ffdatadir','WORK_NACA2408');

NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Generation of a mesh

We use the program NacaFourDigits_Mesh.edp which is a generic mesh generator for classical four-digits NACA profiles of the form NACA_MPXX where M is the camber, P the location of max camber and XX the thickness. This program is designed to be used as follows in shell mode (with M=1, P=4, XX=08):

> FreeFem++ NacaFourDigits_Mesh.edp -M 1 -P 4 -XX 8

Using the driver SF_Mesh, this is interfaced as follows:

mesh = SF_Mesh('NacaFourDigits_Mesh.edp','Options',{'M',1,'P',4,'XX',8});

Plot this mesh :

figure(1);
subplot(2,2,1); SF_Plot(mesh,'title','Mesh (full)');
subplot(2,2,2); SF_Plot(mesh,'title','Mesh (wing)','xlim',[-.5 1],'ylim',[-.75 .75]);
subplot(2,2,3); SF_Plot(mesh,'title','Mesh (Zoom LE)','xlim',[-.4 -.1],'ylim',[-.15 .15]);
subplot(2,2,4); SF_Plot(mesh,'title','Mesh (Zoom TE)','xlim',[.6 .9],'ylim',[-.15 .15]);

In, the sequel we will use the program PotentialFlow.edp which computes the potential flow for given values of incidence $\alpha$ and circulation $\Gamma$ (remind that the potential problem admits multiple solutions parametrized by the circulation $\Gamma$ ).

Compute Potential flow for alpha = 0 ° ; Gamma = 0 (Kutta condition not verified)

alpha = 0; %angle in degrees
Gamma = 0; % circulation
Flow0 = SF_Launch('PotentialFlow.edp','Options',{'alpha',alpha,'Gamma', Gamma},'mesh',mesh,...
                  'DataFile','PotentialFlow.txt','Store','POTENTIALFLOWS');

WARNING -  Folder WORK_NACA2408//POTENTIALFLOWS mentioned for database storing does not exist ; creating it

The Kutta indicator is defined as ux(xc,.001)-ux(xc,-.001) ;

Flow0.Kutta

ans =

   -0.2335

Here this indicator is not zero meaning that the flow is not regular at the trailing edge. See the struture of the flow :

figure(2);

subplot(2,2,[1 2]); SF_Plot(Flow0,'p','xlim',[-.5 1.5],'ylim',[-.75 .75],'colorrange',[-.5 .5],'colormap','redblue','title','Flow for alpha = 0° ; \Gamma = 0');
hold on ; SF_Plot(Flow0,{'ux','uy'},'xlim',[-2 2],'ylim',[-.75 .75]);
subplot(2,2,3); SF_Plot(Flow0,'p','xlim',[-.3 -.2],'ylim',[-.01 .05],'colorrange',[-.5 .5],'cbtitle','p','colormap','redblue','title','(zoom on trailing edge)');
hold on ; SF_Plot(Flow0,{'ux','uy'},'xlim',[-.3 -.2],'ylim',[-.05 .05]);
subplot(2,2,4); SF_Plot(Flow0,'p','xlim',[.7 .8],'ylim',[-.01 .05],'colorrange',[-.5 .5],'cbtitle','p','colormap','redblue','title','(zoom on trailing edge)');
hold on ; SF_Plot(Flow0,{'ux','uy'},'xlim',[.7 .8],'ylim',[-.05 .05]);
pause(0.1);

Computing flow for alpha = 0 degrees with right circulation (Kutta condition verified)

The value of Gamma has been determined by trial-and-error to satisfy the Kutta condition :

Gamma =0.062 ; alpha = 0;
FlowK = SF_Launch('PotentialFlow.edp','Options',{'alpha',alpha,'Gamma', Gamma},'mesh',mesh,...
                  'DataFile','PotentialFlow.txt','Store','POTENTIALFLOWS');

Plot this solution :

figure(4);

subplot(2,2,[1 2]); SF_Plot(FlowK,'p','xlim',[-.5 1.5],'ylim',[-.75 .75],'colorrange',[-.5 .5],'colormap','redblue','title','Flow for alpha = 0° ; \Gamma = \Gamma_K');
hold on ; SF_Plot(Flow0,{'ux','uy'},'xlim',[-2 2],'ylim',[-.75 .75]);
subplot(2,2,3); SF_Plot(FlowK,'p','xlim',[-.3 -.2],'ylim',[-.01 .05],'colorrange',[-.5 .5],'cbtitle','p','colormap','redblue','title','(zoom on trailing edge)');
hold on ; SF_Plot(Flow0,{'ux','uy'},'xlim',[-.3 -.2],'ylim',[-.05 .05]);
subplot(2,2,4); SF_Plot(FlowK,'p','xlim',[.7 .8],'ylim',[-.01 .05],'colorrange',[-.5 .5],'cbtitle','p','colormap','redblue','title','(zoom on trailing edge)');
hold on ; SF_Plot(Flow0,{'ux','uy'},'xlim',[.7 .8],'ylim',[-.05 .05]);
pause(0.1);

This solution has nonzero Lift force, namely

FlowK.Fy

ans =

    0.0620

VISCOUS FLOWS

We consider again alpha = 0 ; we want to ompute a solution for Re = 10^5. We will use the program 'Newton_2D.edp' which computes a steady solution of Navier-Stokes iteration for given parameters Re and alpha.

We first compute a solution for a lower value of Re :

FlowV = SF_Launch('Newton_2D.edp','Options',{'alpha',alpha,'Re', 100},'mesh',mesh,...
                  'DataFile','BaseFlow.txt','Store','VISCOUSFLOWS');

WARNING -  Folder WORK_NACA2408//VISCOUSFLOWS mentioned for database storing does not exist ; creating it

We now progressively increase Re and use continuation method (i.e. the Newton is initiated at each time by the previously computed solution).

At some stages we also use mesh adaptation to optimize the computation.

Re_Tab = [300 600 1000 2000 5000 1e4];
for Re = Re_Tab
  FlowV = SF_Launch('Newton_2D.edp','Options',{'alpha',alpha,'Re', Re},'Init',FlowV,...
                  'DataFile','BaseFlow.txt','Store','VISCOUSFLOWS');
  Mask = SF_Mask(FlowV,[-.3 1 -0.1 0.1 0.005]);
  FlowV = SF_Adapt(FlowV,Mask,'hmin',0.0001,'recompute',false);
end
FlowV = SF_Launch('Newton_2D.edp','Options',{'alpha',alpha,'Re', Re},'Init',FlowV,...
                  'DataFile','BaseFlow.txt','Store','VISCOUSFLOWS') ;

Plot this solution

figure(5) ;hold off; subplot(2,2,[1 2]);
SF_Plot(FlowV,'vort','xlim',[-.5 1.5],'ylim',[-.75 .75],'colorrange',[-10 10],'colormap','redblue');
hold on; SF_Plot(FlowV,{'ux','uy'},'xlim',[-1 3],'ylim',[-.75 .75]);
hold on; SF_Plot(FlowV,'p','contour','only','xlim',[-1 3],'ylim',[-.75 .75],'CStyle','dashedneg');
subplot(2,2,3);
SF_Plot(FlowV,'vort','xlim',[-.4 -.1],'ylim',[-.15 .15],'colorrange',[-10 10],'colormap','redblue');
hold on; SF_Plot(FlowV,{'ux','uy'},'xlim',[-.4 -.1],'ylim',[-.15 .15]);
hold on; SF_Plot(FlowV,'p','contour','only','xlim',[-.4 -.1],'ylim',[-.15 .15],'CStyle','dashedneg');
subplot(2,2,4);
SF_Plot(FlowV,'vort','xlim',[.6 .9],'ylim',[-.15 .15],'colorrange',[-10 10],'colormap','redblue');
hold on; SF_Plot(FlowV,{'ux','uy'},'xlim',[.6 .9],'ylim',[-.15 .15]);
hold on; SF_Plot(FlowV,'p','contour','only','xlim',[.6 .9],'ylim',[-.15 .15],'CStyle','dashedneg');
pause(0.1);

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

figure;
SF_Plot(FlowV,'ux','xlim',[.6 .9],'ylim',[-.15 .15],'colormap','redblue');
hold on; SF_Plot(FlowV,'psi','contour','only','xlim',[.6 .9],'ylim',[-.15 .15],'CStyle','dashedneg');
pause(0.1);

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

Compare the Lift and drag forces in potential and viscous cases

PotentialLift = FlowK.Fy
PotentialDrag = FlowK.Fx
ViscousLift = FlowV.Fy
ViscousDrag = FlowV.Fx

PotentialLift =

    0.0620


PotentialDrag =

   1.0330e-05


ViscousLift =

    0.0127


ViscousDrag =

    0.0166

SUMMARY

SF_Status

% [[PUBLISH]]


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

################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    WORK_NACA2408/
#################################################################################
 
.... CONTENT OF DIRECTORY WORK_NACA2408/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 02:03:15 | 2986    
2     | FFMESH_000002.msh | adapted         | 23-mars-2023 02:04:20 | 8036    
3     | FFMESH_000003.msh | adapted         | 23-mars-2023 02:05:22 | 9498    
4     | FFMESH_000004.msh | adapted         | 23-mars-2023 02:06:39 | 7889    
5     | FFMESH_000005.msh | adapted         | 23-mars-2023 02:07:45 | 10366   
6     | FFMESH_000006.msh | adapted         | 23-mars-2023 02:09:15 | 11907   
7     | FFMESH_000007.msh | adapted         | 23-mars-2023 02:11:01 | 14434   
 
#################################################################################
#################################################################################
 
.... CONTENT OF DIRECTORY WORK_NACA2408/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         | Re         | alpha      | Fx         | Fy        
 1     | FFDATA_000001.txt | MASK         | 23-mars-2023 02:04:17 | FFMESH_000001.msh
 2     | FFDATA_000002.txt | BaseFlow     | 23-mars-2023 02:04:27 | FFMESH_000002.msh | 300        | 0          | 0.10648    | 0.011712  
 3     | FFDATA_000003.txt | MASK         | 23-mars-2023 02:05:19 | FFMESH_000002.msh
 4     | FFDATA_000004.txt | BaseFlow     | 23-mars-2023 02:05:34 | FFMESH_000003.msh | 600        | 0          | 0.072781   | 0.013265  
 5     | FFDATA_000005.txt | MASK         | 23-mars-2023 02:06:34 | FFMESH_000003.msh
 6     | FFDATA_000006.txt | BaseFlow     | 23-mars-2023 02:06:49 | FFMESH_000004.msh | 1000       | 0          | 0.055373   | 0.014166  
 7     | FFDATA_000007.txt | MASK         | 23-mars-2023 02:07:40 | FFMESH_000004.msh
 8     | FFDATA_000008.txt | BaseFlow     | 23-mars-2023 02:07:58 | FFMESH_000005.msh | 2000       | 0          | 0.038379   | 0.015343  
 9     | FFDATA_000009.txt | MASK         | 23-mars-2023 02:09:10 | FFMESH_000005.msh
 10    | FFDATA_000010.txt | BaseFlow     | 23-mars-2023 02:09:32 | FFMESH_000006.msh | 5000       | 0          | 0.023857   | 0.014967  
 11    | FFDATA_000011.txt | MASK         | 23-mars-2023 02:10:55 | FFMESH_000006.msh
 12    | FFDATA_000012.txt | BaseFlow     | 23-mars-2023 02:11:17 | FFMESH_000007.msh | 10000      | 0          | 0.016659   | 0.012357  
 
#################################################################################
 
.... CONTENT OF DIRECTORY WORK_NACA2408/POTENTIALFLOWS
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | alpha      | Gamma      | Kutta      | Fx         | Fy         | Mz        
 1     | FFDATA_000001.txt | PotentialFlow | 23-mars-2023 02:03:24 | FFMESH_000001.msh | 0          | 0          | -0.23349   | -0.00022716 | 0.00018394 | -0.0010077
 2     | FFDATA_000002.txt | PotentialFlow | 23-mars-2023 02:03:35 | FFMESH_000001.msh | 0          | 0.062      | 1.9703e-05 | 1.033e-05  | 0.061957   | 0.013648  
 
#################################################################################
 
.... CONTENT OF DIRECTORY WORK_NACA2408/VISCOUSFLOWS
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | alpha      | Fx         | Fy        
 1     | FFDATA_000001.txt | BaseFlow     | 23-mars-2023 02:03:59 | FFMESH_000001.msh | 100        | 0          | 0.19697    | 0.0094729 
 2     | FFDATA_000002.txt | BaseFlow     | 23-mars-2023 02:04:13 | FFMESH_000001.msh | 300        | 0          | 0.10648    | 0.011712  
 3     | FFDATA_000003.txt | BaseFlow     | 23-mars-2023 02:05:10 | FFMESH_000002.msh | 600        | 0          | 0.072781   | 0.013265  
 4     | FFDATA_000004.txt | BaseFlow     | 23-mars-2023 02:06:24 | FFMESH_000003.msh | 1000       | 0          | 0.055373   | 0.014166  
 5     | FFDATA_000005.txt | BaseFlow     | 23-mars-2023 02:07:31 | FFMESH_000004.msh | 2000       | 0          | 0.038379   | 0.015343  
 6     | FFDATA_000006.txt | BaseFlow     | 23-mars-2023 02:08:57 | FFMESH_000005.msh | 5000       | 0          | 0.023857   | 0.014967  
 7     | FFDATA_000007.txt | BaseFlow     | 23-mars-2023 02:10:37 | FFMESH_000006.msh | 10000      | 0          | 0.016659   | 0.012357  
 8     | FFDATA_000008.txt | BaseFlow     | 23-mars-2023 02:12:13 | FFMESH_000007.msh | 10000      | 0          | 0.016646   | 0.01265   
 
 
#################################################################################

ans = 

  struct with fields:

            MESHES: [1x7 struct]
              MISC: [1x12 struct]
    POTENTIALFLOWS: [1x2 struct]
      VISCOUSFLOWS: [1x8 struct]