FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
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
solver      string 'openfoam'/{'} Use OpenFOAM or default solver
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;
'solver',   'openfoam';
'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.21], 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.p(3,:) = fea.grid.p(3,:) + 0.41;

% Renumber boundaries to match geometry.
fea.grid.b(3,fea.grid.b(3,:)==6) = -7;
fea.grid.b(3,fea.grid.b(3,:)==5) = -8;
fea.grid.b(3,fea.grid.b(3,:)==4) = -9;
fea.grid.b(3,fea.grid.b(3,:)==7) = -10;
[~,ind] = findbdr( fea, 'x<=sqrt(eps)', false );
fea.grid.b(3,ind) = 5;
[~,ind] = findbdr( fea, 'x>=2.5-sqrt(eps)', false );
fea.grid.b(3,ind) = 3;
[~,ind] = findbdr( fea, 'y<=sqrt(eps)', false );
fea.grid.b(3,ind) = 2;
[~,ind] = findbdr( fea, 'y>=0.41-sqrt(eps)', false );
fea.grid.b(3,ind) = 4;
[~,ind] = findbdr( fea, 'z<=sqrt(eps)', false );
fea.grid.b(3,ind) = 1;
[~,ind] = findbdr( fea, 'z>=0.41-sqrt(eps)', false );
fea.grid.b(3,ind) = 6;
fea.grid.b(3,fea.grid.b(3,:)==-7)  = 7;
fea.grid.b(3,fea.grid.b(3,:)==-8)  = 8;
fea.grid.b(3,fea.grid.b(3,:)==-9)  = 9;
fea.grid.b(3,fea.grid.b(3,:)==-10) = 10;
fea.grid.s(:) = 1;
if( strcmp(opt.solver,'openfoam') )
fea.grid = gridrefine( fea.grid, fid );
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 };
fea.phys.ns.sfun            = { opt.sf_u opt.sf_u opt.sf_u opt.sf_p };
if( strcmp(opt.solver,'openfoam') )
[fea.phys.ns.sfun{:}] = deal('sflag1');
end

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

% Parse and solve problem.
fea  = parsephys( fea );
fea  = parseprob( fea );
if( strcmp(opt.solver,'openfoam') )
logfid = fid; if( ~got.fid ), fid = []; end
fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid );
fid = logfid;
else
fea.sol.u = solvestat( fea, 'fid', fid );
end

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

% Error checking.
s_tfx = ['nx*p+',num2str(opt.miu),'*(-2*nx*ux-nz*(uz+vx))'];
s_cd  = ['-2*(',s_tfx,')/(',num2str(opt.rho),'*',num2str(0.2),'^2*',num2str(0.1*0.41),')'];
c_d1  = abs(intbdr(s_cd,fea,[7:10],10));
dp    = evalexpr('p',[0.45 0.55;0.205 0.205;0.21 0.21],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