FEATool Multiphysics  v1.17.0
Finite Element Analysis Toolbox
ex_heattransfer5.m File Reference

Description

EX_HEATTRANSFER5 2D Transient cooling and shrink fitting example.

[ FEA, OUT ] = EX_HEATTRANSFER5( VARARGIN ) This example models transient cooling for shrink fitting of a two part assembly. A tungsten rod chilled to -10 C is inserted into a steel frame heated to 84 C. The assembly is cooled due to convection through a surrounding medium kept at 17 C, and a constant heat transfer coefficient of 50 W/(m^2 K). The time when the maximum temperature has cooled to 70 C should be determined

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.005}         Grid cell size
sfun        string {sflag1}        Finite element shape function
solver      string fenics/{}       Use FEniCS or default solver
iplot       scalar {1}/0           Plot solution (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { 'hmax',     0.005;
             'sfun',     'sflag1';
             'solver',   '';
             'iplot',    1;
             'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});


% Geometry definition.
 r1 = gobj_rectangle( 0, 0.11, 0, 0.12,  'R1' );
 c1 = gobj_circle( [ 0.065 0 ],   0.015, 'C1' );
 c2 = gobj_circle( [ 0.11 0.12 ], 0.035, 'C2' );
 c3 = gobj_circle( [ 0 0.06 ],    0.025, 'C3' );

 r2 = gobj_rectangle( 0.065, 0.16, 0.05, 0.07, 'R2' );
 c4 = gobj_circle( [ 0.065 0.06 ], 0.01, 'C4' );

 fea.geom.objects = { r1 c1 c2 c3 r2 c4 };
 fea = geom_apply_formula( fea, 'R1-C1-C2-C3' );
 fea = geom_apply_formula( fea, 'R2+C4' );


% Grid generation.
 fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', opt.fid );


% Problem definition.
 fea.sdim  = { 'x', 'y' };             % Space coordinate names.
 fea = addphys( fea, @heattransfer );  % Add heat transfer physics mode.
 fea.phys.ht.sfun = { opt.sfun };      % Set shape function.

% Equation coefficients.
 rho_tungsten = 19000;
 cp_tungsten  =   134;
 k_tungsten   =   163;
 rho_steel    =  7500;
 cp_steel     =   470;
 k_steel      =    44;
 fea.phys.ht.eqn.coef{1,end} = { rho_steel rho_tungsten rho_tungsten };   % Density
 fea.phys.ht.eqn.coef{2,end} = { cp_steel  cp_tungsten  cp_tungsten  };   % Heat capcity.
 fea.phys.ht.eqn.coef{3,end} = { k_steel   k_tungsten   k_tungsten   };   % Thermal conductivity.
 fea.phys.ht.eqn.coef{7,end} = { 84 -10 -10 };                            % Initial temperature.

% Boundary conditions.
 n_bdr = max(fea.grid.b(3,:));
 fea.phys.ht.bdr.sel(fea.phys.ht.bdr.sel>0) = 4;
 h_coef = 50;
 for i_bdr=1:n_bdr
   fea.phys.ht.bdr.coef{4,end}{i_bdr}{2} = h_coef;
   fea.phys.ht.bdr.coef{4,end}{i_bdr}{3} = 17;
 end


% Parse physics modes and problem struct.
 fea = parsephys(fea);
 fea = parseprob(fea);


% Compute solution.
 ischeme = 2;
 if( strcmp(opt.solver,'fenics') )
   fea = fenics( fea, 'fid', opt.fid, ...
                 'tstep', 1, 'tmax', 150, 'ischeme', ischeme );
   tlist = fea.sol.t;
 else
   [fea.sol.u, tlist] = solvetime( fea, 'fid', opt.fid, ...
                                   'ischeme', ischeme, ...
                                   'tstep', 1, ...
                                   'tmax', 150, ...
                                   'init', {'T0_ht'} );
 end
 fea.sol.t = tlist;


% Check when max temperature < 70.
 for i=1:length(tlist)
   T_min(i) = min(fea.sol.u(:,i));
   T_max(i) = max(fea.sol.u(:,i));
 end
 ind = find(T_max<70);
 i1 = ind(1);
 i2 = i1 - 1;
 s  = ( T_max(i2) - 70 )/( T_max(i2) - T_max(i1) );
 t_70 = tlist(i2) + s*( tlist(i1) - tlist(i2) );
 u_70 = fea.sol.u(:,i2) + s*( fea.sol.u(:,i1) - fea.sol.u(:,i2) );


% Postprocessing.
 if( opt.iplot>0 )
   figure
   fea_plot = fea;
   fea_plot.sol.u = u_70;
   postplot( fea_plot, 'surfexpr', 'T', 'isoexpr', 'T', 'isolev', 20, 'parent', subplot(1,2,1) )
   title(['Temperature distribution at t = ',num2str(t_70)])

   subplot(1,2,2)
   plot( tlist, T_min, 'b-' )
   hold on
   plot( tlist, T_max, 'r-' )
   grid on
   title('Maximum and minimum temperatures')
   ylabel('Temperature [C]')
   xlabel('Time [s]')
 end


% Error checking.
 out.T_min = T_min;
 out.T_max = T_max;
 out.T_70  = t_70;
 out.pass = abs(out.T_70-138)/138 < 0.02;

 if( nargout==0 )
   clear fea out
 end