FEATool Multiphysics  v1.16.6 Finite Element Analysis Toolbox
ex_navierstokes15.m File Reference

## Description

EX_NAVIERSTOKES15 2D Channel flow driven by periodic pressure difference.

[ FEA, OUT ] = EX_NAVIERSTOKES15( VARARGIN ) Sets up and solves stationary Poiseuille flow in a rectangular channel driven by a periodic pressure difference. The flow profile is constant and should assume a parabolic profile u(y)=dp/dx/2/miu*y*(h-y). 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}           Channel height
l           scalar {2.5}           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.
global dp
dp        = opt.dp;      % Pressure differetial.
% Geometry and grid parameters.
h         = opt.h;       % Height of rectangular domain.
l         = opt.l;       % Length of rectangular domain.
% 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, l, 0, h );
fea.geom.objects = { gobj };
fea.sdim = { 'x' 'y' };   % Coordinate names.

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

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

% Problem definition.
fea = addphys(fea,@navierstokes);        % 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 = parsephys(fea);   % Check and parse physics modes.

% Assign periodic BCs.
[fea.bdr.d{3}{[i_inflow i_outflow]}] = deal( [] );
[fea.bdr.n{3}{[i_inflow i_outflow]}] = deal( @periodic_pressure_bc_2_4 );

% Check and parse problem struct.
fea = parseprob( fea );

% Call to stationary solver.
fea.sol.u = solvestat( fea, 'fid', fid, ...
'nlinasm', [1 1 1], 'tolchg', 1e-2, 'toldef', 1e-4 );

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

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

subplot(3,1,3)
plot( u, p(2,:) )
hold on
plot( u_ref, p(2,:), 'k.' )
legend( 'Computed solution', 'Reference solution','Location','West')
title( 'Velocity profile at x=l/2' )
ylabel( 'y' )
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

%