FEATool Multiphysics  v1.16.5
Finite Element Analysis Toolbox
Shrink Fitting of an Assembly

The following example models shrink fitting of a two part assembly. A chilled tungsten rod is inserted into a slot in a heated steel frame. The rod can be fitted due to thermal expansion of the frame, and when the assembly reaches equilibrium temperature the two parts will be joined due to the ensuing pressure fit.

This model illustrates how to simulate the transient heat transfer process with two domains of different materials. The assembly is cooled due to convection through a surrounding medium kept at Tinf = 17 °C, and the time when the maximum temperature has reached 70 °C should be determined.

FEATool allows for modeling heat transfer through both conduction, that is heat transported by a diffusion process, and also convection, which is heat transported by a fluid due to a flow field. The heat transfer physics mode supports both these processes, and is defined by the following equation

\[ \rho C_p\frac{\partial T}{\partial t} + \nabla\cdot(-k\nabla c) + \rho C_p\mathbf{u}\cdot\nabla T = Q \]

where \(\rho\) is the density, \(C_p\) the heat capacity, \(k\) is the thermal conductivity, \(Q\) heat source term, and \(\mathbf{u}\) a vector valued convective velocity field.

Here, the surrounding cooling medium is not modeled directly, the convective term in the heat transfer equation can therefore be dropped, and the cooling effects are incorporated into the model by using natural convection boundary conditions instead. By not directly modeling the surrounding medium significant computational savings can be made.

Tutorial

This section describes how to set up and solve the thermal shrink fitting example with the FEATool graphical user interface (GUI). The model is available as an automated tutorial by selecting Model Examples and Tutorials... > Heat Transfer > Shrink Fitting of an Assembly from the File menu, viewed as a videotutorial", or alternatively, follow the step-by-step instructions below. Note that the CFDTool interface differ slightly from the FEATool Multiphysics instructions described in the following.

  1. To start a new model click the New Model toolbar button, or select New Model... from the File menu.
  2. Select the Heat Transfer physics mode from the Select Physics drop-down menu. (Note that for CFDTool the physics selection is done in the Equation settings dialog box.)

  3. Press OK to finish the physics mode selection.

The model geometry consists of a cross section of the frame and rod. The frame can be generated by first creating a rectangle, and then removing three circles and the space for the slot.

  1. To create a rectangle, first click on the Create square/rectangle Toolbar button. Then left click in the main plot axes window, and hold down the mouse button. Move the mouse pointer to draw the shape outline, and release the button to finalize the shape.
  2. Select R1 in the geometry object Selection list box.
  3. To modify and edit the selected rectangle, click on the Inspect/edit selected geometry object Toolbar button to open the Edit Geometry Object dialog box.
  4. Enter 0 into the xmin edit field.
  5. Enter 0.11 into the xmax edit field.
  6. Enter 0 into the ymin edit field.
  7. Enter 0.12 into the ymax edit field.

  8. Press OK to finish and close the dialog box.

Now create the three circles which are to be removed from R1.

  1. To create a circle or ellipse, first click on the Create circle/ellipse Toolbar button. Then left click in the main plot axes window, and hold down the mouse button. Move the mouse pointer to draw the shape outline, and release the button to finalize the shape.
  2. Select E1 in the geometry object Selection list box.
  3. To modify and edit the selected ellipse, click on the Inspect/edit selected geometry object Toolbar button to open the Edit Geometry Object dialog box.
  4. Enter 0.065 0 into the center edit field.
  5. Enter 0.015 into the xradius edit field.
  6. Enter 0.015 into the yradius edit field.

  7. Press OK to finish and close the dialog box.
  8. To create a circle or ellipse, first click on the Create circle/ellipse Toolbar button. Then left click in the main plot axes window, and hold down the mouse button. Move the mouse pointer to draw the shape outline, and release the button to finalize the shape.
  9. Select E2 in the geometry object Selection list box.
  10. To modify and edit the selected ellipse, click on the Inspect/edit selected geometry object Toolbar button to open the Edit Geometry Object dialog box.
  11. Enter 0.11 0.12 into the center edit field.
  12. Enter 0.035 into the xradius edit field.
  13. Enter 0.035 into the yradius edit field.

  14. Press OK to finish and close the dialog box.
  15. To create a circle or ellipse, first click on the Create circle/ellipse Toolbar button. Then left click in the main plot axes window, and hold down the mouse button. Move the mouse pointer to draw the shape outline, and release the button to finalize the shape.
  16. Select E3 in the geometry object Selection list box.
  17. To modify and edit the selected ellipse, click on the Inspect/edit selected geometry object Toolbar button to open the Edit Geometry Object dialog box.
  18. Enter 0 0.06 into the center edit field.
  19. Enter 0.025 into the xradius edit field.
  20. Enter 0.025 into the yradius edit field.

  21. Press OK to finish and close the dialog box.

  22. Select Combine Objects... from the Geometry menu.
  23. Enter R1 - E1 - E2 - E3 into the Geometry Formula edit field.

  24. Press OK to finish and close the dialog box.

Now we need to create the geometry for the rod by joining a smaller rectangle and a circle.

  1. To create a rectangle, first click on the Create square/rectangle Toolbar button. Then left click in the main plot axes window, and hold down the mouse button. Move the mouse pointer to draw the shape outline, and release the button to finalize the shape.
  2. Select R2 in the geometry object Selection list box.
  3. To modify and edit the selected rectangle, click on the Inspect/edit selected geometry object Toolbar button to open the Edit Geometry Object dialog box.
  4. Enter 0.065 into the xmin edit field.
  5. Enter 0.16 into the xmax edit field.
  6. Enter 0.05 into the ymin edit field.
  7. Enter 0.07 into the ymax edit field.

  8. Press OK to finish and close the dialog box.
  9. To create a circle or ellipse, first click on the Create circle/ellipse Toolbar button. Then left click in the main plot axes window, and hold down the mouse button. Move the mouse pointer to draw the shape outline, and release the button to finalize the shape.
  10. Select E4 in the geometry object Selection list box.
  11. To modify and edit the selected ellipse, click on the Inspect/edit selected geometry object Toolbar button to open the Edit Geometry Object dialog box.
  12. Enter 0.065 0.06 into the center edit field.
  13. Enter 0.01 into the xradius edit field.
  14. Enter 0.01 into the yradius edit field.

  15. Press OK to finish and close the dialog box.
  16. Select R2 and E4 in the geometry object Selection list box.

  17. Press the + / Add geometry objects Toolbar button.

To make a slot in the frame for the rod, first copy it and subtract the copy from the frame.

  1. Select CJ1 in the geometry object Selection list box.
  2. Press the Copy and/or transform selected geometry object Toolbar button.
  3. Select the Make copy of geometry object check box.

  4. Press OK to finish and close the dialog box.
  5. Select CS1 and CJ1 in the geometry object Selection list box.

  6. Press the - / Subtract geometry objects Toolbar button.

The final geometry consists of two geometry objects, one domain for the frame, and one for the rod.

  1. Switch to Grid mode by clicking on the corresponding Mode Toolbar button.

The default grid may be too coarse to ensure an accurate solution. Decreasing the grid size and generating a finer grid can resolve curved boundaries better.

  1. Enter 0.0025 into the Grid Size edit field and press the Generate button to call the grid generation algorithm.

  2. Switch to Equation mode by clicking on the corresponding Mode Toolbar button.

Equation and material coefficients are specified in Equation/Subdomain mode. In the Equation Settings dialog box enter the coefficients for density, heat capacity, thermal conductivity, and initial temperature for each material. (Note that the Equation Settings dialog box may look different for CFDTool.)

  1. Select 1 in the Subdomains list box.
  2. Enter rho_tungsten into the Density edit field.
  3. Enter cp_tungsten into the Heat capacity edit field.
  4. Enter k_tungsten into the Thermal conductivity edit field.
  5. Enter -10 into the Initial condition for T edit field.

  6. Select 2 in the Subdomains list box.
  7. Enter rho_steel into the Density edit field.
  8. Enter cp_steel into the Heat capacity edit field.
  9. Enter k_steel into the Thermal conductivity edit field.
  10. Enter 84 into the Initial condition for T edit field.

  11. Press OK to finish the equation and subdomain settings specification.

The Model Constants and Expressions functionality can be used to define and store convenient expressions, which then are available in the point, equation, boundary coefficients, and as postprocessing expressions. Here it is used to define the material parameters.

  1. Press the Constants Toolbar button, or select the corresponding entry from the Equation menu, and enter the following variables in the Model Constants and Expressions dialog box. Press Enter after the last expression or use the Add Row button to expand the expression list.
Name Expression
rho_tungsten 19000
cp_tungsten 134
k_tungsten 163
rho_steel 7500
cp_steel 470
k_steel 44
h_coef 50


  1. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.
  2. In the Boundary Settings dialog box, select the Heat flux boundary condition for all the boundaries. Enter h_coef in the edit field for the convective transfer coefficient h, and also enter 17 for the surrounding reference temperature Tinf.

  3. Now that the problem is fully specified, press the Solve mode button to switch to solve mode. Since this is a time dependent study, open the solver settings and select the Time-Dependent solver. Set the Time step to 1, Simulation time to 150, and then press Solve to start the solution process.

After the solver has finished, the toolbox will automatically switch to postprocessing mode, and show the resulting temperature field. Add some contour levels to more clearly see the different temperature levels.

  1. Press the Plot Options Toolbar button.
  2. Select the Contour Plot check box.
  3. Enter 20 into the Number or specified vector of contour levels to plot edit field.

  4. Press Apply to plot and visualize the selected postprocessing options.

Plot solutions at different times and verify that the assembly has cooled to a temperature of 70 degrees around the time t = 138 s. Note that both the colorbar and limits field (in the Postprocessing dialog box) will show the minimum and maximum surface plot values.

  1. Select 138 from the Available solutions/times drop-down menu, and press OK to plot and visualize the temperature at that time.

More advanced processing can also be done by exporting the simulation data to the command line interface, where it can be accessed, modified, and used in scripts. (See for example the section on CLI postprocessing below.)

The shrink fitting of an assembly heat transfer model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB m-script text file (available as the example ex_heattransfer5 script file), or GUI script (.fes) file.

CLI Postprocessing

More advanced postprocessing can be done on the MATLAB command line interface (CLI) console. To export the model data to the MATLAB CLI, select Export Model Data Struct to MATLAB option from the File menu. The model data will now be available as a variable named fea. The following script commands will find the time and plot the solution when the maximum temperature has reached 70 degrees.

% Find minimum, maxium, and time 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) );

% Plot temperature field at t_70.
fig = 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)])

% Plot min/max temperature curves as a function of time.
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]')

Video