![]() | FEATool
v1.7.1 Finite Element Analysis Toolbox |
EX_NAVIERSTOKES6 2D incompressible time-dependent flow around a cylinder.
[ FEA, OUT ] = EX_NAVIERSTOKES6( VARARGIN ) Benchmark example for time-dependent flow around a cylinder.
[1] Schaefer M, Turek S. Benchmark computations of laminar flow around a cylinder. In Flow Simulation with High-Performance Computers II, Notes on Numerical Fluid Dynamics, vol. 52. Vieweg, 1996; 547-566.
[2] 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
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 %