|
FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
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
iplot scalar 0/{1} Plot solution and error (=1)
.
Output Value/(Size) Description
-----------------------------------------------------------------------------------
fea struct Problem definition struct
out struct Output struct
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, 'ddtSchemes', 'CrankNicolson 0.9', 'deltaT', opt.dt, 'endTime', opt.tmax, 'nproc', 1, 'fid', fid, 'logfid', logfid );
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