FEATool Multiphysics
v1.17.1
Finite Element Analysis Toolbox
|
EX_NATURAL_CONVECTION 2D Example for natural convection of air in a square cavity.
[ FEA, OUT ] = EX_NATURAL_CONVECTION( VARARGIN ) Sets up and solves a natural convection benchmark problem. Reference solutions are for example reported in G. Davis "Natural convection of air in a square cavity a bench mark numerical solution", IJNMF vol. 3, 249-254 (1983).
Accepts the following property/value pairs.
Input Value/{Default} Description ----------------------------------------------------------------------------------- Ra scalar {1e3} Rayleigh number Pr scalar {0.71} Prandtl number l scalar {1} Side length of cavity igrid scalar 0/{1} Cell type (0=quadrilaterals, 1=triangles) hmax scalar {0.05} Max grid cell size sf_u string {sflag2} Shape function for velocity sf_p string {sflag1} Shape function for pressure sf_T string {sflag2} Shape function for temperature iphys scalar 0/{1} Use physics mode to define problem (=1) solver string {} Solver selection default, openfoam, fenics iplot scalar 0/{1} Plot solution and error (=1) . Output Value/(Size) Description ----------------------------------------------------------------------------------- fea struct Problem definition struct out struct Output struct
cOptDef = { ... 'Ra', 1e3; 'Pr', 0.71; 'l', 1; 'igrid', 0; 'hmax', 0.05; 'sf_u', 'sflag2'; 'sf_p', 'sflag1'; 'sf_T', 'sflag2'; 'iphys', 1; 'solver', ''; 'iplot', 1; 'tol', 0.1; 'fid', 1 }; [got,opt] = parseopt(cOptDef,varargin{:}); fid = opt.fid; % Reference values. Ra_ref = [ 1e3 1e4 1e5 1e6 ]; u_max_ref = [ 3.649 16.178 34.73 64.63 ]; y_max_ref = [ 0.813 0.823 0.855 0.850 ]; v_max_ref = [ 3.697 19.617 68.59 219.36 ]; x_max_ref = [ 0.178 0.119 0.066 0.0379 ]; Nu_mean_ref = [ 1.118 2.243 4.519 8.800 ]; % Model parameters. Ra = opt.Ra; % Rayleigh number. Pr = opt.Pr; % Prandtl number. l = opt.l; % Length of rectangular domain. sf_u = opt.sf_u; % FEM shape function type for velocity. sf_p = opt.sf_p; % FEM shape function type for pressure. sf_T = opt.sf_T; % FEM shape function type for temperature. % Geometry definition. fea.sdim = { 'x' 'y' }; fea.geom.objects = { gobj_rectangle( 0, l, 0, l ) }; % Grid generation. if( opt.igrid<=0 ) fea.grid = rectgrid( round(l/opt.hmax), round(l/opt.hmax), [0 l;0 l] ); if( opt.igrid<0 ) fea.grid = quad2tri(fea.grid); end else fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid ); end % Boundary conditions. n_bdr = max(fea.grid.b(3,:)); % Increment number of boundaries. dtol = opt.hmax; ib_l = findbdr( fea, ['x<',num2str(dtol)] ); % Right boundary number. ib_r = findbdr( fea, ['x>',num2str(l-dtol)] ); % Left boundary number. ib_b = findbdr( fea, ['y<',num2str(dtol)] ); % Bottom boundary number. ib_t = findbdr( fea, ['y>',num2str(l-dtol)] ); % Top boundary number. % Add pressure point constraint on point closest to origin. [~,ix] = min( fea.grid.p(1,:).^2 + fea.grid.p(2,:).^2 ); fea.pnt.index = ix; fea.pnt.type = 'constr'; fea.pnt.dvar = 'p'; fea.pnt.expr = 0'; % Problem definition. if ( opt.iphys==1 || got.solver ) fea = addphys(fea,@navierstokes); % Add Navier-Stokes equations physics mode. fea.phys.ns.eqn.coef{2,end} = { Pr }; fea.phys.ns.eqn.coef{4,end} = { [num2str(Ra*Pr),'*T'] }; if( any(strcmp(opt.solver,{'openfoam'})) ) fea.phys.ns.sfun = { 'sflag1', 'sflag1', 'sflag1' }; else fea.phys.ns.sfun = { sf_u sf_u sf_p }; % Set shape functions. end fea = addphys(fea,@heattransfer); % Add heat transfer physics mode. if( any(strcmp(opt.solver,{'openfoam'})) ) fea.phys.ht.sfun = { 'sflag1' }; else fea.phys.ht.sfun = { sf_T }; end fea.phys.ht.eqn.coef{4,end} = { fea.phys.ns.dvar{1} }; fea.phys.ht.eqn.coef{5,end} = { fea.phys.ns.dvar{2} }; fea.phys.ht.bdr.sel([ib_l ib_r]) = 1; fea.phys.ht.bdr.coef{1,end}{ib_l} = 1; fea = parsephys(fea); % Check and parse physics modes. else fea.dvar = { 'u' 'v' 'p' 'T' }; % Dependent variable name. fea.sfun = { sf_u sf_u sf_p sf_T }; % Shape function. % Define equation system. cvelx = fea.dvar{1}; % Convection velocities. cvely = fea.dvar{2}; fea.eqn.a.form = { [2 3 2 3;2 3 1 1] [2;3] [1;2] []; [3;2] [2 3 2 3;2 3 1 1] [1;3] []; [2;1] [3;1] [] []; [] [] [] [2 3 2 3;2 3 1 1] }; fea.eqn.a.coef = { {2*Pr Pr cvelx cvely} Pr -1 []; Pr {Pr 2*Pr cvelx cvely} -1 []; 1 1 [] []; [] [] [] {1 1 cvelx cvely} }; fea.eqn.f.form = { 1 1 1 1 }; fea.eqn.f.coef = { 0 [num2str(Ra*Pr),'*',fea.dvar{4}] 0 0 }; % Define boundary conditions. fea.bdr.d = cell(4,n_bdr); [fea.bdr.d{1:2,:}] = deal(0); fea.bdr.d{4,ib_l} = 1; fea.bdr.d{4,ib_r} = 0; fea.bdr.n = cell(4,n_bdr); [fea.bdr.n{4,[ib_t ib_b n_bdr]}] = deal(0); [fea.bdr.n{3,setdiff(1:n_bdr,n_bdr)}] = deal(0); [fea.bdr.n{1:2,n_bdr}] = deal(0); end % Parse and solve problem. fea = parseprob(fea); % Check and parse problem struct. if( strcmp(opt.solver,'openfoam') ) logfid = fid; if( ~got.fid ), fid = []; end fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid, 'nproc', 1 ); fid = logfid; elseif( strcmp(opt.solver,'fenics') ) fea = fenics( fea, 'fid', fid ); else fea.sol.u = solvestat(fea,'fid',fid); % Call to stationary solver. end % Postprocessing. if ( opt.iplot>0 ) figure subplot(1,2,1) postplot(fea,'surfexpr','sqrt(u^2+v^2)') title('Velocity field') subplot(1,2,2) postplot(fea,'surfexpr','T') title('Temperature') end % Error checking. out.err = nan; out.pass = nan; iref = find( Ra==Ra_ref ); if ( ~isempty(iref) ) n_evalution_points = 3*l/opt.hmax; x_eval = linspace( 0, l, n_evalution_points ); x_mid = l/2*ones( 1, n_evalution_points ); u_eval = evalexpr( 'u', [x_mid; x_eval], fea ); [u_max,ix] = max( u_eval ); y_max = x_eval( ix ); v_eval = evalexpr( 'v', [x_eval; x_mid], fea ); [v_max,ix] = max( v_eval ); x_max = x_eval( ix ); Nu_mean = abs( intbdr( 'Tx', fea, 1, 2 ) + intbdr( 'Tx', fea, 2, 2 ) )/2; out.u_max = u_max; out.y_max = y_max; out.v_max = v_max; out.x_max = x_max; out.Nu_mean = Nu_mean; out.err = [ abs(u_max_ref(iref)-u_max)/u_max_ref(iref); abs(y_max_ref(iref)-y_max)/y_max_ref(iref); abs(v_max_ref(iref)-u_max)/v_max_ref(iref); abs(x_max_ref(iref)-x_max)/x_max_ref(iref); abs(Nu_mean_ref(iref)-Nu_mean)/Nu_mean_ref(iref) ]; out.pass = all( out.err<opt.tol ); end if ( nargout==0 ) clear fea out end