FEATool  v1.7.1
Finite Element Analysis Toolbox
 All Files Functions Pages
ex_navierstokes13.m File Reference

Description

EX_NAVIERSTOKES13 3D Example flow past a cylinder.

[ FEA, OUT ] = EX_NAVIERSTOKES13( VARARGIN ) Sets up and solves stationary and laminar 3D flow past a cylinder. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
rho         scalar {1}             Density
miu         scalar {0.001}         Molecular/dynamic viscosity
uin         scalar {0.45}          Magnitude of inlet velocity
sf_u        string {sf_hex_Q1nc}   Shape function for velocity
sf_p        string {sf_disc0}      Shape function for pressure
iplot       scalar 0/{1}           Plot solution and error (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = {   ...
   'rho',      1;
   'miu',      0.001;
   'uin',      0.45;
   'nlev',     1;
   'sf_u',     'sf_hex_Q1nc';
   'sf_p',     'sf_disc0';
   'tol',      0.2;
   'iplot',    1;
   'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


% Geometry.
 fea.sdim = { 'x', 'y', 'z' };
 gobj1 = gobj_block( 0, 2.5, 0, 0.41, 0, 0.41, 'B1' );
 gobj2 = gobj_cylinder( [0.5 0 0.2], 0.05, 0.41, 2, 'C1' );
 fea.geom.objects = { gobj1 gobj2 };
 fea = geom_apply_formula( fea, 'B1-C1' );


% Grid generation.
 ns = 8*2^(opt.nlev-1);
 r  = [0.05 0.06 0.08 0.11 0.15];
 x  = [0.41 0.5 0.7 1 1.4 1.8 2.2] + 0.3;
 for ilev=2:opt.nlev
   r = sort( [ r (r(1:end-1)+r(2:end))/2 ] );
   x = sort( [ x (x(1:end-1)+x(2:end))/2 ] );
 end

 grid1 = ringgrid( r, 4*ns, [], [], [0.5;0.2] );
 grid2 = holegrid( ns, 2^(opt.nlev-1), [0.3 0.71;0 0.41], 0.15, [0.5;0.2] );
 grid2 = gridmerge( grid1, 5:8, grid2, 1:4 );
 grid3 = rectgrid( x, ns, [0.71 2.5;0 0.41] );
 fea.grid = gridmerge( grid3, 4, grid2, 6 );
 grid4 = rectgrid( 1, ns, [0 0.3;0 0.41] );
 fea.grid = gridmerge( fea.grid, 10, grid4, 2 );
 fea.grid = gridextrude( fea.grid, 5, 0.41, 2 );
 fea.grid = gridscale( fea.grid, [1 1 -1] );
 fea.grid = sort_bdr( fea.grid, fea.geom, [] );
 [~,ix] = findbdr( fea , 'z<=sqrt(eps)', 0 );
 fea.grid.b(3,ix) = 1;
 [~,ix] = findbdr( fea , 'z>=(0.41-sqrt(eps))', 0 );
 fea.grid.b(3,ix) = 6;


% 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_u opt.sf_p };


% Boundary conditions.
 fea.phys.ns.bdr.sel(5) = 2;
 fea.phys.ns.bdr.sel(3) = 3;
 fea.phys.ns.bdr.coef{2,end}{1,5} = ['16*',num2str(opt.uin),'*(y*z*(0.41-y)*(0.41-z))/0.41^4'];


% Parse and solve problem.
 fea  = parsephys( fea );
 fea  = parseprob( fea );
 data = featflow_data( 3 );
 data.nlmax = 2;
 fea  = featflow( fea, 'mode', 'export', 'data', data, 'fid', fid );
 fea  = featflow( fea, 'mode', 'solve',  'fid', fid );
 pause(3)
 fea  = featflow( fea, 'mode', 'import', 'fid', fid );


% Postprocessing.
 if( opt.iplot>0 )
   postplot( fea, 'sliceexpr', 'sqrt(u^2+v^2+w^2)' )
 end


% Error checking.
 s_tfx = ['nx*p+',num2str(miu),'*(-2*nx*ux-nz*(uz+vx))'];
 s_cd  = ['-2*(',s_tfx,')/(',num2str(rho),'*',num2str(0.2),'^2*',num2str(0.1*0.41),')'];
 c_d1  = intbdr(s_cd,fea,[7:10],10);
 dp    = evalexpr('p',[0.35 0.55;0.2 0.2;0.2 0.2],fea);
 out.err = [abs(c_d1-6.185333)/6.185333 ...
            abs(dp(1)-dp(2)-0.171007)/0.171007];
 out.pass = all(out.err < [0.15 0.6]);


 if ( nargout==0 )
   clear fea out
 end