FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
ex_navierstokes8.m File Reference

## Description

EX_NAVIERSTOKES8 2D Example for axisymmetric incompressible stationary flow in a constricted circular pipe.

[ FEA, OUT ] = EX_NAVIERSTOKES8( VARARGIN ) Sets up and solves stationary axisymmetric Poiseuille flow in a constricted circular pipe. The inflow profile is constant and the outflow should assume a parabolic profile. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
rho         scalar {1}             Density
miu         scalar {1}             Molecular/dynamic viscosity
uin         scalar {1}             Inflow velocity (constant/mean)
l           scalar {3}             Channel length
hmax        scalar {0.1}           Max grid cell size
sf_u        string {sflag1}        Shape function for velocity
sf_p        string {sflag1}        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

ex_navierstokes8b

# Code listing

 cOptDef = { ...
'rho',      1;
'miu',      1;
'uin',      1;
'r',        1;
'l',        3;
'hmax',     0.05;
'sf_u',     'sflag1';
'sf_p',     'sflag1';
'solver',   '';
'iplot',    1;
'fid',      1 };
[got,opt] = parseopt( cOptDef, varargin{:} );
fid       = opt.fid;

% Geometry definition.
r = opt.r;   % Pipe radius.
l = opt.l;   % Pipe length.
gobj1 = gobj_rectangle( 0, r,   0,     l*2/3, 'R1' );
gobj2 = gobj_rectangle( 0, r/2, l*2/3, l,     'R2' );
gobj3 = gobj_circle( [r l*2/3], r/2, 'C1' );
fea.geom.objects = { gobj1 gobj2 gobj3 };
fea = geom_apply_formula( fea, 'R1+R2-C1' );
fea.sdim = { 'r' 'z' };   % Space coordinate names.

% Grid generation.
fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );

% Boundary specifications.
i_inflow   = 1;       % Inflow boundary number.
i_outflow  = 5;       % Outflow boundary number.
i_symmetry = [3 6];   % Symmetry boundary numbers.

% Problem definition.
fea = addphys( fea, {@navierstokes 1} );   % Add Navier-Stokes equations physics mode.
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 };

% Boundary conditions.
dtol = 1e-3;
i_in  = findbdr( fea, ['z<=',num2str(dtol)] );
i_out = findbdr( fea, ['z>=',num2str(3-dtol)] );
i_sym = findbdr( fea, ['r<=',num2str(dtol)] );
fea.phys.ns.bdr.sel(i_in)   = 2;
fea.phys.ns.bdr.sel(i_out)  = 3;
fea.phys.ns.bdr.coef{2,end}{2,i_in} = opt.uin;
fea.phys.ns.bdr.sel(i_sym) = 5;
fea = parsephys(fea);  % Parse physics mode.

% Parse and solve problem.
fea = parseprob(fea);   % Check and parse problem struct.
if( strcmp(opt.solver,'openfoam') )
logfid = fid; if( ~got.fid ), fid = []; end
fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid );
fid = logfid;
else
jac.form  = {[1;1] [1;1] [];[1;1] [1;1] []; [] [] []};
jac.coef  = {'r*rho_ns*ur' 'r*rho_ns*uz' []; 'r*rho_ns*wr' 'r*rho_ns*wz' []; [] [] []};
fea.sol.u = solvestat( fea, 'fid', fid, 'nsolve', 2, 'jac', jac );
end

% Error checking.
r = linspace( 0, opt.r/2, 20 );
z = 0.9*opt.l*ones( 1, 20 );
U = evalexpr( 'sqrt(u^2+w^2)', [r;z], fea )';
u_fac = 4;   % Due to contraction to 1/2 radius.
U_ref = 2*opt.uin*u_fac*( 1 - ( r/(opt.r/2) ).^2 );
err = sqrt( sum((U-U_ref).^2)/sum(U_ref.^2) );

% Postprocessing.
if( opt.iplot>0 )
figure
subplot(1,3,1)
postplot( fea, 'surfexpr', 'sqrt(u^2+w^2)', 'arrowexpr', {'u' 'w'} )
hold on
plot( r, z, 'k--' )
title( 'Velocity field' )

subplot(1,3,2)
postplot( fea, 'surfexpr', 'p', 'evaltype', 'exact', 'isoexpr', 'p' )
title( 'Pressure' )

subplot(1,3,3)
plot( r, U,     'b-' )
hold on
plot( r, U_ref, 'r-' )
title('Velcity profile at z=0.9*l')
end