Description

EX_NAVIERSTOKES14 Axisymmetric flow in a pipe due to pressure difference.

[ FEA, OUT ] = EX_NAVIERSTOKES14( VARARGIN ) Sets up and solves stationary Poiseuille flow in an axisymmetric pipe driven by a pressure difference. The flow profile is constant and should assume a parabolic profile u(r)=dp/dx/2/miu*(r+h/2)*(h/2-r). Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
rho         scalar {0.1}           Density
miu         scalar {0.2}           Molecular/dynamic viscosity
dp          scalar {0.3}           Pressure difference between in and out-flow
h           scalar {0.5}           Pipe/channel height/width
l           scalar {2.5}           Pipe/channel length
igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
hmax        scalar {0.04}          Max grid cell size
sf_u        string {sflag2}        Shape function for velocity
sf_p        string {sflag1}        Shape function for pressure
iphys       scalar 0/{1}           Use physics mode to define problem (=1)
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',      0.1;
             'miu',      0.2;
             'dp',       0.3;
             'h',        0.5;
             'l',        2.5;
             'igrid',    1;
             'hmax',     0.5/10;
             'sf_u',     'sflag2';
             'sf_p',     'sflag1';
             'iphys',    1;
             'iplot',    1;
             'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


% Model parameters.
 rho       = opt.rho;     % Density.
 miu       = opt.miu;     % Molecular/dynamic viscosity.
 dp        = opt.dp;      % Pressure differetial.
% Geometry and grid parameters.
 h         = opt.h;       % Height/width of pipe.
 l         = opt.l;       % Length of pipe.
% Discretization parameters.
 sf_u      = opt.sf_u;    % FEM shape function type for velocity.
 sf_p      = opt.sf_p;    % FEM shape function type for pressure.


% Geometry definition.
 gobj = gobj_rectangle( 0, h/2, 0, l );
 fea.geom.objects = { gobj };
 fea.sdim = { 'r' 'z' };   % Coordinate names.


% Grid generation.
 fea.grid = rectgrid(round(l/opt.hmax),round(h/2/opt.hmax),[0 h/2;0 l]);
 if( opt.igrid~=0 )
   fea.grid = quad2tri( fea.grid );
 end


% Boundary conditions.
 dtol      = opt.hmax/2;
 i_inflow  = findbdr( fea, ['z<',num2str(dtol)] );     % Inflow boundary number.
 i_outflow = findbdr( fea, ['z>',num2str(l-dtol)] );   % Outflow boundary number.
 i_axis    = findbdr( fea, ['r<',num2str(dtol)] );     % Symmetry axis boundary number.


% Problem definition.
 fea = addphys(fea,@navierstokes,'ns',true);   % Add Navier-Stokes equations physics mode.
 fea.phys.ns.eqn.coef{1,end} = { rho };
 fea.phys.ns.eqn.coef{2,end} = { miu };
 fea.phys.ns.sfun = { sf_u sf_u sf_p };   % Set shape functions.


 fea.phys.ns.bdr.sel([i_inflow i_outflow]) = 4;
 fea.phys.ns.bdr.coef{4,end}{3,i_inflow}   = dp;         % Set inflow pressure.
 fea.phys.ns.bdr.coef{4,end}{3,i_outflow}  = 0;          % Set outflow pressure.
 fea.phys.ns.bdr.sel(i_axis)               = 5;

 fea       = parsephys(fea);                 % Check and parse physics modes.
 fea       = parseprob( fea );               % Check and parse problem struct.
 fea.sol.u = solvestat( fea, 'fid', fid );   % Call to stationary solver.


% Postprocessing.
 s_velm = 'sqrt(u^2+v^2)';
 s_refsol = [num2str(dp),'/',num2str(l),'/2/',num2str(miu),'*(r+',num2str(h),'/2)*(',num2str(h),'/2-r)'];   % Definition of velocity profile.
 p = [ linspace(0,h/2,25); l/2*ones(1,25) ];
 u = evalexpr( s_velm, p, fea );
 u_ref = evalexpr( s_refsol, p, fea );
 if ( opt.iplot>0 )
   figure
   subplot(1,3,1)
   postplot(fea,'surfexpr',s_velm,'evaltype','exact')
   title('Velocity field')

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

   subplot(1,3,3)
   plot( u, p(1,:) )
   hold on
   plot( u_ref, p(1,:), 'k.' )
   legend( 'Computed solution', 'Reference solution','Location','West')
   title( 'Velocity profile at z=l/2' )
   ylabel( 'r' )
 end


% Error checking.
 err = sqrt(sum((u-u_ref).^2)/sum(u_ref.^2));
 if( ~isempty(fid) )
   fprintf(fid,'\nL2 Error: %f\n',err)
   fprintf(fid,'\n\n')
 end


 out.err  = err;
 out.pass = err<0.05;
 if ( nargout==0 )
   clear fea out
 end


%