FEATool Multiphysics  v1.16.6 Finite Element Analysis Toolbox
ex_heattransfer10.m File Reference

## Description

EX_HEATTRANSFER10 Conjugate heat transfer test example with multiple domains.

[ FEA, OUT ] = EX_HEATTRANSFER10( VARARGIN ) Sets up and solves a multiple domain conjugate heat transfer example.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.0005}        Max grid cell size
sf_u        string {sflag2}        Shape function for velocity
sf_p        string {sflag1}        Shape function for pressure
sf_T        string {sflag1}        Shape function for temperature
solver      string openfoam/fenics/{} Use OpenFOAM, FEniCS, 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 = { 'hmax',     1/10;
'sf_u',     'sflag1';
'sf_p',     'sflag1';
'sf_T',     'sflag1';
'solver',   '';
'iplot',    1;
'tol',      0.1;
'fid',      1 };
[got,opt] = parseopt( cOptDef, varargin{:} );
fid       = opt.fid;

% Geometry definition.
fea.sdim         = { 'x' 'y' };
fea.geom.objects = { gobj_rectangle(0, 1, 0, 0.5), ...
gobj_rectangle(0, 1, 0.5, 1, 'R2') };

% Grid generation.
n = round(1/opt.hmax);
fea.grid = gridmerge( rectgrid(2*n, n, [0, 1; 0, 0.5]), 3, rectgrid(2*n, n, [0, 1; 0.5, 1]), 1);
fea.grid.s( selcells(fea,'y >= 0.5') ) = 2;
fea.grid = gridbdrx(fea.grid);

% Problem definition.
fea.phys.ns.eqn.coef{1,end} = {1, 1000};
fea.phys.ns.eqn.coef{2,end} = {1, 959e-6};
fea.phys.ns.prop.active(:,1) = 0;
fea.phys.ns.bdr.sel(6) = 2;
fea.phys.ns.bdr.coef{2,end}{1,6} = 0.00001;
fea.phys.ns.bdr.sel(4) = 4;

fea.phys.ht.eqn.coef{1,end} = {8000, 1000};
fea.phys.ht.eqn.coef{2,end} = {450, 4181};
fea.phys.ht.eqn.coef{3,end} = {80, 4181*959e-6/6.62};
fea.phys.ht.eqn.coef{4,end} = {0, 'u'};
fea.phys.ht.eqn.coef{5,end} = {0, 'v'};
fea.phys.ht.bdr.sel([1,6]) = 1;
fea.phys.ht.bdr.coef{1,end}{1} = 400;
fea.phys.ht.bdr.coef{1,end}{6} = 300;
fea.phys.ht.bdr.sel([2,3,5]) = 3;
fea.phys.ht.eqn.coef{end}{1} = 300;
fea.phys.ht.eqn.coef{end}{2} = 300;

% Parse and solve problem.
fea = parsephys(fea);
fea = parseprob(fea);

if( strcmp(opt.solver,'fenics') )
fea = fenics( fea, 'fid', fid );
elseif( strcmp(opt.solver,'openfoam') )
% fea.sol.u = openfoam( fea, 'fid', fid );
fea.sol.u = openfoam( fea, 'fid', fid, 'ddtScheme', 'CrankNicolson', 'endTime', 1e5, 'maxDeltaT', 1e5/100 );
else
fea.sol.u = solvestat( fea, 'fid', fid, 'maxnit', 50 );
end

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

subplot(1,2,2)
postplot( fea, 'surfexpr', 'T' )
title('Temperature')
end

% Average temperature at outlet.
out.ref = [1.4111e-5, 309.7998, 394.3849];
out.val = [evalexpr('u',[1;0.75],fea), ...
evalexpr('T',[1;0.75],fea), ...
evalexpr('T',[1;0.5],fea)];
out.err = abs(out.val - out.ref)./out.ref;
out.pass = all(out.err < opt.tol);

if( nargout==0 )
clear fea out
end