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

## Description

EX_NAVIERSTOKES5 2D Vortex flow with analytical solution.

[ FEA, OUT ] = EX_NAVIERSTOKES5( VARARGIN ) Sets up and solves time dependent flow of a decaying vortex for which an analytical solution is known. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
nu          scalar {0.01}          Kinematic viscosity
k           scalar {pi/2}          Vortex number
xdim        {[-1 3]}               Domain min/max x-coordinates
ydim        {[-1 3]}               Domain min/max y-coordinates
igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
hmax        scalar {0.2}           Max grid cell size
dt          scalar {0.01}          Time step size
tmax        icalar {0.1}           Simluation duration
ischeme     scalar {3}             Time stepping scheme
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


# Code listing

 cOptDef = { ...
'nu',       0.01;
'k',        pi/2;
'xdim',     [-1 3];
'ydim',     [-1 3];
'igrid',    0;
'hmax',     0.2;
'dt',       0.01;
'tmax',     0.1;
'ischeme',  2;
'sf_u',     'sflag1';
'sf_p',     'sflag1';
'solver',   '';
'iplot',    1;
'tol',      [0.03, 0.03, 0.12];
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Exact solutions.
k = num2str(opt.k);
u_ref = ['-cos(',k,'*x).*sin(',k,'*y)*exp(-',num2str(2*opt.nu*opt.k^2),'*t)'];
v_ref = [' sin(',k,'*x).*cos(',k,'*y)*exp(-',num2str(2*opt.nu*opt.k^2),'*t)'];
p_ref = ['-1/4*(cos(2*',k,'*x)+cos(2*',k,'*y))*exp(-4*',num2str(opt.nu*opt.k^2),'*t)'];

% Grid generation.
fea.sdim = { 'x' 'y' };
if ( opt.igrid==0 )
fea.grid = rectgrid( round(abs(diff(opt.xdim))/opt.hmax), round(diff(abs(opt.ydim))/opt.hmax), [opt.xdim;opt.ydim] );
else
fea.geom.objects = { gobj_rectangle( opt.xdim(1), opt.xdim(2), opt.ydim(1), opt.ydim(2) ) };
fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );
end

% Problem definition.
fea = addphys( fea, @navierstokes );
fea.phys.ns.eqn.coef{1,end} = { 1 };
fea.phys.ns.eqn.coef{2,end} = { opt.nu };
fea.phys.ns.eqn.coef{5,end} = { u_ref };
fea.phys.ns.eqn.coef{6,end} = { v_ref };
fea.phys.ns.eqn.coef{7,end} = { p_ref };
fea.phys.ns.bdr.sel = [2 2 2 2];
[fea.phys.ns.bdr.coef{2,end}{1,:}] = deal( u_ref );
[fea.phys.ns.bdr.coef{2,end}{2,:}] = deal( v_ref );
[fea.phys.ns.bdr.coef{2,end}{3,:}] = deal( p_ref );
fea.phys.ns.sfun = { opt.sf_u opt.sf_u opt.sf_p };
fea = parsephys( fea );

% Remove integratl constraints.
if( isfield(fea,'constr') )
fea = rmfield(fea,'constr');
end

% Parse and solve problem.
fea = parseprob( fea );
[fea.bdr.d{1}{:}]  = deal( u_ref );   % Exact boundary conditions.
[fea.bdr.d{2}{:}]  = deal( v_ref );
[fea.bdr.d{3}{:}]  = deal( p_ref );
fea.pnt = struct;
fea.bdr.n = cell(3,max(fea.grid.b(3,:)));
if( strcmp(opt.solver,'openfoam') )
logfid = fid; if( ~got.fid ), fid = []; end
[fea.sol.u,fea.sol.t] = openfoam( fea, 'fid', fid, 'logfid', logfid, 'ddtScheme', 'CrankNicolson', 'deltaT', opt.dt, 'endTime', opt.tmax );
fid = logfid;
else
[fea.sol.u,fea.sol.t] = solvetime( fea, 'fid', fid, 'init', { u_ref v_ref p_ref }, 'ischeme', opt.ischeme, 'tstep', opt.dt, 'tmax', opt.tmax, 'nstbwe', 0, 'icub', 2 );
end

% Postprocessing.
if ( opt.iplot>0 )
subplot(1,2,1)
postplot( fea, 'surfexpr', 'sqrt(u^2+v^2)', 'arrowexpr', {'u', 'v'}, 'arrowcolor', 'k' )
title('Velocity field')

subplot(1,2,2)
postplot( fea, 'surfexpr', 'p' )
title('Pressure')
end

% Error checking.
x = fea.grid.p(1,:)';
y = fea.grid.p(2,:)';
n = length(fea.sol.t);
for i_sol=1:n
t = fea.sol.t(i_sol);

u_i = evalexpr( 'u', fea.grid.p, fea, i_sol );
v_i = evalexpr( 'v', fea.grid.p, fea, i_sol );
p_i = evalexpr( 'p', fea.grid.p, fea, i_sol );

u_r = eval( u_ref );
v_r = eval( v_ref );
p_r = eval( p_ref );

errnm(i_sol,1) = norm( u_i - u_r )/norm( u_r );
errnm(i_sol,2) = norm( v_i - v_r )/norm( v_r );
errnm(i_sol,3) = norm( p_i - p_r )/norm( p_r );
end
out.err  = errnm;
out.pass = all( errnm(:) < reshape(repmat(opt.tol,n,1),[],1) );

if ( nargout==0 )
clear fea out
end