FEATool Multiphysics  v1.13
Finite Element Analysis Toolbox
ex_navierstokes6.m File Reference

Description

EX_NAVIERSTOKES6 2D incompressible time-dependent flow around a cylinder.

[ FEA, OUT ] = EX_NAVIERSTOKES6( VARARGIN ) Benchmark example for time-dependent flow around a cylinder.

References

[1] John V. Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder. International Journal for Numerical Methods in Fluids 2004; 44:777-788.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
rho         scalar {1}             Density
miu         scalar {0.001}         Molecular/dynamic viscosity
igrid       scalar {2}             Grid type: >0 regular (igrid refinements)
                                              <0 unstruc. grid (with hmax=|igrid|)
dt          scalar {0.02}          Time step size
instatbc    boolean {true}         Use instationary boundary conditions
ischeme     scalar {2}             Time stepping scheme
solver      string {}              Solver selection default, openfoam, su2
sf_u        string {sflag2}        Shape function for velocity
sf_p        string {sf_disc1}      Shape function for pressure
iplot       scalar 0/{1}           Plot solution (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { ...
   'rho',      1;
   'miu',      0.001;
   'igrid',    2;
   'dt',       0.02;
   'instatbc'  true;
   'ischeme'   2;
   'solver',   '';
   'sf_u',     'sflag2';
   'sf_p',     'sf_disc1';
   'iplot',    1;
   'tol',      0.1;
   'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;
 if( any(strcmp(opt.solver,{'openfoam','su2'})) )
   opt.instatbc = false;
 end


% Grid generation.
 fea.sdim = { 'x' 'y' };
 if( opt.igrid>=1 )
   fea.grid = cylbenchgrid( opt.igrid );
   fea.grid.s(:) = 1;
 else
   gobj1 = gobj_rectangle( 0, 2.2, 0, 0.41, 'R1' );
   gobj2 = gobj_circle( [0.2 0.2], 0.05, 'C1' );
   fea.geom.objects = { gobj1 gobj2 };
   fea = geom_apply_formula( fea, 'R1-C1' );
   fea.grid = gridgen( fea, 'hmax', abs(opt.igrid), 'fid', fid );
 end


% Boundary conditions.
 inflow_boundary  = findbdr( fea, ['x<=',num2str(sqrt(eps))] );   % Inflow boundary number.
 outflow_boudnary = findbdr( fea, ['x>=',num2str(2.1)] );         % Outflow boundary number.
 if( opt.instatbc )
   inflow_bc = '6*sin(pi*t/8)*(y*(0.41-y))/0.41^2';
 else
   inflow_bc = '6*(y*(0.41-y))/0.41^2';
 end


% Problem definition.
 fea = addphys( fea, @navierstokes );
 fea.phys.ns.eqn.coef{1,end} = { opt.rho };
 fea.phys.ns.eqn.coef{2,end} = { opt.miu };
 if( any(strcmp(opt.solver,{'openfoam','su2'})) )
   fea.phys.ns.sfun = { 'sflag1', 'sflag1', 'sflag1' };
 else
   fea.phys.ns.sfun = { opt.sf_u, opt.sf_u, opt.sf_p };
 end
 fea.phys.ns.bdr.sel( inflow_boundary  ) = 2;
 fea.phys.ns.bdr.sel( outflow_boudnary ) = 4;
 fea.phys.ns.bdr.coef{2,end}{1,inflow_boundary} = inflow_bc;
 fea = parsephys( fea );


% Call solver
 if( strcmp(opt.solver,'openfoam') )
   if( opt.ischeme==1 )
     ddtScheme = 'backward';
   else
     ddtScheme = 'CrankNicolson 0.9';
   end
   logfid = fid; if( ~got.fid ), fid = []; end
   [fea.sol.u,tlist] = openfoam( fea, 'ddtScheme', ddtScheme, 'deltaT', opt.dt, 'endTime', 8, 'fid', fid, 'logfid', logfid );
   fid = logfid;
 elseif( strcmp(opt.solver,'su2') )
   logfid = fid; if( ~got.fid ), fid = []; end
   [fea.sol.u,tlist] = su2( fea, 'fid', fid, 'logfid', logfid, 'tstep', opt.dt, 'tmax', 8, 'ischeme', opt.ischeme, 'wrtfreq', floor(8/opt.dt/400) );
   fid = logfid;
 else
   [fea.sol.u,tlist] = solvetime( fea, 'fid', fid, 'tstep', opt.dt, 'tmax', 8, 'maxnit', 5, 'ischeme', opt.ischeme );
 end

% Benchmark quantities.
 [ c_d, c_l, dp ] = calc_bench_quant( fea, 1 );
 [c_d_max, i] = max( c_d );
 t_c_d_max    = tlist( i );
 [c_l_max, i] = max( c_l );
 t_c_l_max    = tlist( i );
 out.c_d = c_d;
 out.c_l = c_l;
 out.dp  = dp;
 out.c_d_max = c_d_max;
 out.c_l_max = c_l_max;
 out.t_c_d_max = t_c_d_max;
 out.t_c_l_max = t_c_l_max;
 if( opt.instatbc )
   [~, i] = min(abs(tlist-8));
   dp_t8  = dp(i);
   out.dp_t8 = dp_t8;

   comp_data = [ t_c_d_max c_d_max t_c_l_max c_l_max dp_t8 ];
   ref_data  = [ 3.93625 2.950923849 5.6925 0.47834818 -0.11162153 ];
   out.err   = abs(ref_data-comp_data)./abs(ref_data);
   out.pass  = all( out.err < opt.tol );
 else
   found = false;
   for j=i+1:length(c_l)
     if( j==length(c_l) )
       break
     end
     if( c_l(j) > c_l(j-1) && c_l(j) > c_l(j+1) )
       found = true;
       c_l_max2 = c_l(j);
       t_c_l_max2 = tlist(j);
       break
     end
   end
   for j=i-1:-1:1
     if( c_l(j) > c_l(j-1) && c_l(j) > c_l(j+1) )
       found = true;
       c_l_max2 = c_l(j);
       t_c_l_max2 = tlist(j);
       break
     end
   end
   f = 1/(t_c_l_max2-t_c_l_max);
   St = 0.1*f/1.5;

   out.St = St;
   t_dp = t_c_l_max + 1/2/f;
   [~,j] = find( tlist > t_dp, 1 );
   dp_f2 = dp(j-1) + (t_dp - tlist(j-1))/(tlist(j) - tlist(j-1)) * (dp(j)-dp(j-1));
   out.dp_f2 = dp_f2;

   c_d(tlist<0.5) = 0;
   [c_d_max, i] = max( c_d );
   comp_data = [ c_d_max c_l_max St dp_f2 ];
   ref_data  = [ 3.23 1.00 0.3 2.48 ];
   out.err   = abs(ref_data-comp_data)./abs(ref_data);
   out.pass  = all( out.err([1,2,4]) < opt.tol );
 end


 if( ~isempty(fid) )
   fmtf = '%12.8f |';
   fmts = '%12s |';
   fmt  = ['|      ',repmat(fmtf,1,4+double(opt.instatbc)),'\n'];
   fmts = ['|      ',repmat(fmts,1,4+double(opt.instatbc)),'\n'];
   fmtl = ['|------',repmat('-------------+',1,4+double(opt.instatbc))];
   fmtl = [fmtl(1:end-1),'|\n'];
   fprintf( fid, '\n\n' );
   fprintf( fid, fmtl );
   if( opt.instatbc )
     fprintf( fid, fmts, 't(cd_max)', 'cd_max', 't(cl_max)', 'cl_max', 'dp(t=8)' );
     fprintf( fid, fmtl );
     fprintf( fid, fmt, t_c_d_max, c_d_max, t_c_l_max, c_l_max, dp_t8 );
   else
     fprintf( fid, fmts, 'cd_max', 'cl_max', 'St', 'dp' );
     fprintf( fid, fmtl );
     fprintf( fid, fmt, c_d_max, c_l_max, St, dp_f2 );
   end
   fprintf( fid, fmtl );
   fmt = ['| Ref. ',repmat(fmtf,1,4+double(opt.instatbc)),'\n'];
   fprintf( fid, fmt, ref_data );
   fprintf( fid, fmtl );
   fprintf( fid, '\n\n' );
 end


% Postprocessing.
 if( opt.iplot>0 )
   figure
   fea = parseprob( fea );

   subplot(2,2,1)
   postplot( fea, 'surfexpr', 'sqrt(u^2+v^2)', 'arrowexpr', {'u' 'v'}, 'arrowcolor', 'k' )
   title('Velocity field at t=8')

   ix = tlist > 1;
   if( any(ix) )
     tlist = tlist(ix);
     c_d = c_d(ix);
     c_l = c_l(ix);
     dp  = dp(ix);
   end

   subplot(2,2,3)
   plot( tlist, c_d )
   title('drag coefficient')

   subplot(2,2,4)
   plot( tlist, c_l )
   title('lift coefficient')

   subplot(2,2,2)
   plot( tlist, dp )
   title('pressure difference')
 end


 if( nargout==0 )
   clear fea out
 end


%