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