 |
FEATool Multiphysics
v1.16.5
Finite Element Analysis Toolbox
|
|
FEATool Multiphysics
v1.16.5
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)
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; ...
'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 )
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'] };
fea.phys.ns.sfun = { sf_u sf_u sf_p }; % Set shape functions.
fea = addphys(fea,@heattransfer); % Add heat transfer physics mode.
fea.phys.ht.sfun = { sf_T };
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.
fea.sol.u = solvestat(fea,'fid',fid); % Call to stationary solver.
% 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