|
FEATool Multiphysics
v1.17.5
Finite Element Analysis Toolbox
|
|
|
FEATool Multiphysics
v1.17.5
Finite Element Analysis Toolbox
|
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
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, 'nproc', 1 );
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