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
ischeme     scalar {3}             Time stepping scheme
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; ...
   'ischeme'   3; ...
   'sf_u',     'sflag2'; ...
   'sf_p',     'sf_disc1'; ...
   'iplot',    1; ...
   'itest',    0; ...
   'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


% Grid generation.
 fea.sdim = { 'x' 'y' };
 if( opt.igrid>=1 )
   fea.grid = cylbenchgrid( opt.igrid );
 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.
 inflow_bc = '6*sin(pi*t/8)*(y*(0.41-y))/0.41^2';


% Problem definition.
 fea = addphys( fea, @navierstokes );
 fea.phys.ns.eqn.coef{1,end} = { opt.rho };
 fea.phys.ns.eqn.coef{2,end} = { opt.miu };
 fea.phys.ns.sfun = { opt.sf_u opt.sf_u opt.sf_p };
 fea.phys.ns.bdr.sel( inflow_boundary  ) = 2;
 fea.phys.ns.bdr.sel( outflow_boudnary ) = 3;
 fea.phys.ns.bdr.coef{2,end}{1,inflow_boundary} = inflow_bc;
 fea = parsephys( fea );


% Call solver
 [fea.sol.u,tlist] = solvetime( fea, 'fid', fid, 'tstep', opt.dt, 'tmax', 8*(~opt.itest), 'maxnit', 1+4*(~opt.itest), 'ischeme', opt.ischeme );


% 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 );
 [~, i]       = min(abs(tlist-8));
 dp_t8        = dp(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;
 out.dp_t8 = dp_t8;

 if( ~isempty(fid) )
   fmtf = '%12.8f |';
   fmts = '%12s |';
   fmt  = ['|      ',repmat(fmtf,1,5),'\n'];
   fmts = ['|      ',repmat(fmts,1,5),'\n'];
   fmtl = ['|------',repmat('-------------+',1,5)];
   fmtl = [fmtl(1:end-1),'|\n'];
   fprintf( fid, '\n\n' );
   fprintf( fid, fmtl );
   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 );
   fprintf( fid, fmtl );
   ref_data = [ 3.93625 2.950923849 5.6925 0.47834818 -0.11162153 ];
   fmt  = ['| Ref. ',repmat(fmtf,1,5),'\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')

   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


%