Finite Element Analysis Toolbox
Multi-Simulation Heat Exchanger

This model simulates how the flow of cool airflow is heated while moving through a tube-fin heat exchanger. Due to several symmetry planes only a small section of the heat exchanger geometry actually needs to be simulated, as illustrated in the following image.

This model illustrates the multi-solver simulation process for a coupled tube and fin heat exchanger. The flow field is first solved for using the OpenFOAM CFD solver, after which the temperature field is computed with the built-in FEA solver (using the known flow field as a input to the heat equation).

The process splits and separates the equations, effectively turning a one-way coupled model to a two step solution process with two sequential problems. This allows for significant savings in computational time and resources, since the best and most efficient physics solvers can be used for each individual sub-problem.

Tutorial

This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Multiphysics > Multi-Simulation Heat Exchanger from the File menu, viewed as a video tutorial, or alternatively, follow the step-by-step instructions below. (Note that the external OpenFOAM solver needs to be installed for this model to run correctly.)

  1. To start a new model click the New Model toolbar button, or select New Model... from the File menu.
  2. Select the 3D radio button.
  3. Select the Navier-Stokes Equations physics mode from the Select Physics drop-down menu.

  4. Press OK to finish the physics mode selection.

A heat transfer physics mode for the temperature field will be added and coupled later, after the flow field has been computed.

By utilizing the symmetry in the y and z directions the computational domain of the airflow can be significantly reduced to a slice between two fins and one tube. The resulting geometry can be constructed by subtracting the fins and a cylinder from a block. The first step is to create the main block for the interior of the domain.

  1. Press the Create cube/block Toolbar button.
  2. Enter 20 into the xmax edit field.
  3. Enter 5 into the ymax edit field.
  4. Press OK to finish and close the dialog box.

Then create a cylinder and subtract it from the block.

  1. Press the Create cylinder/cone Toolbar button.
  2. Enter 10 0 0 into the edit field for the center of the cylinder.
  3. Enter 2.5 into the edit fields for both radius1 and radius2.
  4. Enter 0 0 1 into the axis edit field.
  5. Press OK to finish and close the dialog box.
  6. Select B1 and C1 in the geometry object Selection list box.
  7. Press the - / Subtract geometry objects Toolbar button.

Create the lower fin, and then make a copy with a z-translation to move it to the upper side.

  1. Press the Create cube/block Toolbar button.
  2. Enter 5 into the xmin edit field.
  3. Enter 15 into the xmax edit field.
  4. Enter 5 into the ymax edit field.
  5. Enter 0.0625 into the zmax edit field.
  6. Press OK to finish and close the dialog box.
  7. Select B2 in the geometry object Selection list box.
  8. Press the Copy/transform selected geometry object Toolbar button.
  9. Enter 1 into the Number of copies to make edit field.
  10. Enter 0 0 0.9375 into the Displacement vector (x, y, and z-components) edit field.
  11. Press OK to finish and close the dialog box.

Lastly, remove the two fins using the geometry formula CS1 - B2 - B3.

  1. Select Combine Objects... from the Geometry menu.
  2. Enter CS1-B2-B3 into the Geometry Formula edit field.

  3. Press OK to finish and close the dialog box. The completed geometry should then look like the following.

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

Generate a grid with the maximum target mesh size set to 0.2. Although this is a rather coarse mesh, it saves computational time and is good enough for demonstration purposes and as a first study.

  1. Enter 0.2 into the Grid Size edit field.

  2. Press the Generate button to call the grid generation algorithm.

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

Enter a non-dimensionalized unit density of 1, and a viscosity of 0.00526. This is equivalent to a Reynolds number of 190 based on the distance between the heat exchanger fins.

  1. Enter 1 into the Density edit field.

  2. Enter 0.00526 into the Viscosity edit field.

  3. Press OK to finish the equation and subdomain settings specification.
  4. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.

First, set the velocity on the inlet boundary in the x-direction to 1.

  1. Select 5 in the Boundaries list box.
  2. Select Inlet/velocity from the Navier-Stokes Equations drop-down menu.
  3. Enter 1 into the Velocity in x-direction edit field.

Then select the Outflow/pressure condition for the outlet boundary.

  1. Select 11 in the Boundaries list box.
  2. Select Outflow/pressure from the Navier-Stokes Equations drop-down menu.

  3. Select the Wall/no-slip condition for the boundaries representing the fins and the cylinder (boundaries 3, 7-10, 13, 14, and 16).

Finally, select the Symmetry/slip condition for the rest of the boundaries.

  1. Select 1, 2, 4, 6, 12, 15, and 17 in the Boundaries list box.

  2. Select Symmetry/slip from the Navier-Stokes Equations drop-down menu.

  3. Press OK to finish the boundary condition specification.
  4. Switch to Solve mode by clicking on the corresponding Mode Toolbar button.

The OpenFOAM CFD solver will be used to first solve for the flow field. Open the OpenFOAM solver settings dialog box and set the tolerance for convergence to 1e-4.

  1. Press the OpenFOAM Toolbar button.

  2. Enter 1e-4 into the Stopping criteria/tolerance for initial residuals edit field.

  3. Press the Solve button to start the OpenFOAM solver. The view will switch to show the convergence process for the solution variables.

After the flow problem has been solved FEATool will automatically switch to postprocessing mode and display the resulting velocity field.

Open the Postprocessing settings dialog box and change from surface to slice plot to help see the interior of the flow field.

  1. Press the Plot Options Toolbar button.
  2. Clear the Enable/disable surface plot check box.
  3. Select the Enable/disable slice plot check box.

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

One can now clearly see that there is a large wake behind the cylinder, and how the fins create a very thin low velocity boundary layer.

To couple and study the temperature field, switch back to Equation mode and add a Heat Transfer physics mode to the model.

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

  2. First deactivate the equation for the flow field by de-selecting the active button. This means that the flow variables will not be solved for and held constant, which saves computational effort. (Note that this decoupling is only possible for one-way coupled multiphysics problems. If the flow field and properties also depend on the temperature, both physics modes must be fully coupled and solved for together.)

  3. Next switch to the tab with a + symbol, where additional physics modes can be added.
  4. Select the Heat Transfer physics mode from the Select Physics drop-down menu.
  5. Press the Add Physics >>> button.

Set the non-dimensionalized thermal conductivity to 3.76e-3, while leaving the density and heat capacity at their default unit values. This is equivalent to a Prandtl number of 1.4 (indicating a slight weighting toward convective transport).

  1. Enter 3.76e-3 into the Thermal conductivity edit field.

To couple the flow field to the convective terms for the temperature, enter the dependent variable names u, v, and w in the corresponding edit fields.

  1. Enter u into the Convection velocity in x-direction edit field.
  2. Enter v into the Convection velocity in y-direction edit field.

  3. Enter w into the Convection velocity in z-direction edit field.

As this model is dominated by convective flow effects some amount of artificial and numerical stabilization is appropriate to add, which will ease convergence and smooth out oscillations.

  1. Press the Artificial Stabilization button.
  2. Select the Enable/disable streamline diffusion check box.

  3. Enter 1 into the Streamline diffusion tuning parameter edit field.

  4. Press OK to finish and close the dialog box.
  5. Press OK to finish the equation and subdomain settings specification.
  6. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.

For the temperature boundary conditions set the inlet temperature to 0, and the surfaces of the surrounding fins and cylinder to 1.

  1. Switch to the ht tab.
  2. Select 5 in the Boundaries list box.
  3. Select Temperature from the Heat Transfer drop-down menu.
  4. Enter 0 into the Temperature edit field.
  5. Select 3, 7-10, 13, 14, and 16 in the Boundaries list box.
  6. Select Temperature from the Heat Transfer drop-down menu.
  7. Enter 1 into the Temperature edit field.

For the outflow boundary select Convective flux/outflow.

  1. Select 11 in the Boundaries list box.
  2. Select Convective flux/outflow from the Heat Transfer drop-down menu.

And finally select Thermal insulation/symmetry for the symmetry boundaries.

  1. Select 1, 2, 4, 6, 12, 15, and 17 in the Boundaries list box.

  2. Select Thermal insulation/symmetry from the Heat Transfer drop-down menu.

  3. Press OK to finish the boundary condition specification.
  4. Switch to Solve mode by clicking on the corresponding Mode Toolbar button.

Press the Restart button to solve the problem for the temperature field using the already computed flow field (as the Navier-Stokes equations physics mode was deactivated earlier). (Do not use the usual = solve button, as this would clear the already computed flow field and instead re-compute the initial conditions.)

  1. Press the Restart/Solve with last computed solution as initial condition Toolbar button to start the solution process for the temperature.

After the solution process is done the temperature field can now be plotted and visualized.

  1. Press the Plot Options Toolbar button.
  2. Select Temperature, T from the Predefined slice plot expressions drop-down menu.

  3. Press OK to plot and visualize the selected postprocessing options.

One can clearly see how the fluid is heated by both the cylinder and walls, and transported away in the direction of the flow.

The multi-simulation heat exchanger multiphysics 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_heat_exchanger2 script file), or GUI script (.fes) file.

CLI Postprocessing

To visualize the full 3D solution the model data can be exported to the MATLAB command line interface (CLI) console with the Export Model Data Struct to MATLAB option from the File menu. The following advanced script can then be used to recreate the full model 

% Create and plot full heat exchanger geometry.
geom.objects = {};
zoffset = 0;
for i=1:21
  tag = ['B',num2str(i)];
  gobj = gobj_block( [5], [15], [0], [40], [zoffset-0.0625], [zoffset+0.0625], tag );
  geom.objects{i} = gobj;
  zoffset = zoffset + 1;
end

yoffset = 0;
for i=1:4
  tag = ['C',num2str(i)];
  gobj = gobj_cylinder( [10 yoffset+5 -5], [2.5], [30], [0 0 1], tag );
  geom.objects = [ geom.objects, {gobj} ];
  yoffset = yoffset + 10;
end

fig = figure;
hold on
plotgeom( geom, 'parent', fig, 'labels', 'off' )

% Plot simulation data.
fea.grid.p(3,:) = fea.grid.p(3,:) + 10;
fea.grid.p(2,:) = fea.grid.p(2,:) + 5;
h = postplot(fea, 'parent', fig, 'sliceexpr', 'T', 'slicex', [], 'slicey', [], 'slicez', 10+0.3);
hp = findall(h,'type','patch');
delete(setdiff(h,hp))

pdata = get( hp, {'vertices', 'faces', 'facevertexcdata'} );
p = pdata{1};
for yoff=[10, 10, 10]
  p(:,2) = p(:,2) + yoff;
  patch( 'vertices', p, 'faces', pdata{2}, 'facevertexcdata', pdata{3},  'facecolor', 'interp', 'linestyle', 'none' );
end

p = pdata{1};
p(:,2) = 10 - p(:,2);
for yoff=[0, 10, 10, 10]
  p(:,2) = p(:,2) + yoff;
  patch( 'vertices', p, 'faces', pdata{2}, 'facevertexcdata', pdata{3},  'facecolor', 'interp', 'linestyle', 'none' );
end

% Offset fea data.
fea.grid.p(2,:) = fea.grid.p(2,:) + 10;
hfea = postplot( fea, 'parent', fig, 'streamexpr', {'u', 'v', 'w'} );

% Camera settings.
view(3)
axis off
alpha(0.05)
for i=length(geom.objects)-3:length(geom.objects)
  set(findall(0,'tag',['fea_geomp',num2str(i)]),'FaceAlpha',0.2)
end
set(findall(0,'tag','geom_slice'),'FaceAlpha',1)
set(hp,'FaceAlpha',1)
delete(findall(0,'type','light'))
camlight headlight

Video