The following post shows how to define and solve a two parameter thermal FEA design and optimization model problem with FEATool Multiphysics and MATLAB’s optimization functionality.

Design Optimization for Cooling a Heated Block with Target Temperature

In addition to featuring an easy-to-use graphical user interface (GUI), FEATool Multiphysics also includes native application programming interface (API) for FEA and physics simulation fully compatible with MATLAB functions and toolboxes such as Simulink and Optimization. Simulation models can first be prototyped and defined in the GUI, and later exported to equivalent stand alone simulation scripts. These in turn can easily be modified to run parametric simulation studies, or integrated into more advanced applications, such as here for design optimization.

In the following, FEATool’s heat transfer solver is combined with MATLAB’s fminbnd optimization function to solve a practical thermal design problem, that is finding the optimal number and simultaneously the diameter of water-filled cooling pipes embedded in a heated block, so that the peak temperature reaches a prescribed target value. The approach presented here requires no additional toolboxes beyond a base MATLAB installation.

Problem Description - Geometry and Mesh

The geometry consists of a heated block (with dimensions 400 × 200 × 20 mm) mounted on top of a cooling plate with the same footprint and thickness. Between 1 and 9 cylindrical cooling pipes of variable diameter (5-19 mm) run lengthwise through the cooling plate, carrying chilled water at a mass flow rate of 0.02 kg/s. A constant volumetric heat source of 500 W is applied to the top block, and natural convection to ambient air (20 W/(m²·K), 298.15 K) acts on all external surfaces. The goal is to find the pipe count n and diameter D that minimize |T − T_max| subject to the constraint T_max ≤ T, where the target temperature is T^* = 358.15 K (85°C).

geom_2D = l_block_cooling_geometry( n_pipes, diameter, 2 );

grid_2D = gridgen( geom_2D, 'hmax', 0.0025 );
grid_3D = gridextrude( grid_2D, 20, 0.4, -2 );

By exploiting the cross-sectional symmetry of the configuration, only one half of the domain needs to be modeled. The 2D cross-sectional geometry, comprising the top block, the cooling plate, and up to [n/2] circular pipe cut-outs, is generated with the FEATool CAD geometry API meshed using the built-in mesh generation function gridgen with a maximum element size of 2.5 mm. The resulting 2D mesh is then extruded 20 cells lengths along the 400 mm length to produce a full 3D tetrahedral mesh.

Physics and Boundary Conditions

The heat transfer physics mode is used to simulate thermal conduction and convection. Material properties and source terms are assigned per subdomain, the top block has a thermal conductivity of 80 W/(m·K) and includes the volumetric heat source, while the bottom cooling plate (subdomain 2) has thermal conductivity 0.5 W/(m·K).

All external faces are assigned a natural convection boundary condition using the heat flux boundary condition type (4) with a heat transfer coefficient of 20 W/(m²·K) and surrounding ambient temperature 298.15 K.

fea.phys.ht.bdr.sel(external_boundaries) = 4;
[fea.phys.ht.bdr.coef{4,end}{external_boundaries}] = deal({0, 20, 298.15, 0, 0});

For the pipe surfaces, the forced convection heat transfer coefficient is computed from the flow conditions using the Dittus-Boelter correlation (Nusselt number Nu) when the pipe Reynolds number is turbulent (Re ≥ 10 000), or the constant laminar value of Nu = 3.66 otherwise. And the cooler water temperature is assumed to be 278.15 K.

Re = (4 * mass_flow_rate) / (pi * diameter * viscosity);
if( Re >= 10000 )
  Nu = 0.023 * Re^0.8 * Pr^0.4;
else
  Nu = 3.66;
end
k_thermal_conductivity = 0.598;
c_ht = (Nu * k_thermal_conductivity) / diameter;

fea.phys.ht.bdr.sel(pipe_boundaries) = 4;
[fea.phys.ht.bdr.coef{4,end}{pipe_boundaries}] = deal({0, c_ht, 278.15, 0, 0});

The symmetry planes are assigned boundary condition type 3 (Thermal insulation/symmetry with zero heat flux).

fea.phys.ht.bdr.sel(symmetry_boundaries) = 3;

Lastly, the physics mode and FEA model struct is parsed, and solved with the stationary built-in multiphysics FEA solver solvestat The maximum temperature can be obtained using the minmaxsubd function.

fea = parsephys( fea );
fea = parseprob( fea );

fea.sol.u = solvestat( fea );

[~, Tmax] = minmaxsubd( 'T', fea );

Optimization Strategy

This optimization problem has one integer variable (n, the pipe count) and one continuous variable (D, the pipe diameter). Rather than treating this as a single black-box mixed-integer problem the structure is exploited by decomposing it into an outer loop over all integer values of n, and an inner scalar minimization over D for each fixed n using fminbnd. This nested approach is both efficient and requires no additional toolboxes.

The objective passed to fminbnd is a penalized function that combines the primary objective |T^* − T_max| with a large penalty weight (1e4) for constraint violation:

function [ cost ] = penObj( D, n, Tstar, penaltyWeight, cache )

ff   = cachedFEA( n, D, Tstar, cache );
cost = ff.Fval + penaltyWeight * max(0, ff.Ineq);

To avoid redundant FEA solver calls, which dominate the computational cost, results are stored and cached for later retrieval. Any repeated evaluation of the same (n, D) pair immediately returns the cached result:

function [ ff ] = cachedFEA( n, D, Tstar, cache )

key = sprintf('%d_%.8f', n, D);
if( isKey(cache, key) )
  ff = cache(key);

else
  ff = l_run_fea_optimization_problem( [n; D], Tstar );
  cache(key) = ff;
end

A warm-start strategy further reduces the number of function evaluations, the optimal diameter found for the previous pipe count, n, is used as the center of a narrowed search bracket for the next pipe count n+1, since adjacent pipe counts tend to require similar pipe diameters. A fallback to the full diameter range is triggered if the warm-started bracket fails to find a feasible solution.

bracketHalf = (Dmax - Dmin) * 0.3;
a = max(Dmin, D0 - bracketHalf);
b = min(Dmax, D0 + bracketHalf);
[Dopt, cost] = fminbnd( f, a, b, opts );

% Fallback to full range if constraint is violated.
if( ff.Ineq > 1e-2 )
  [Dopt2, cost2] = fminbnd( f, Dmin, Dmax, opts );
  if( cost2 < cost )
    Dopt = Dopt2;  cost = cost2;
  end
end

Results

Running the script sweeps all pipe counts from 1 to 9 and reports the optimal diameter, maximum temperature, penalized cost, and total cumulative FEA evaluations for each, typical results are shown below.

n    D (m)      Tmax (K)     Cost         FEA calls
---------------------------------------------------
1    0.0190     421.64       634961.9589  23
2    0.0190     408.12       499702.6643  44
3    0.0190     386.45       283014.0704  65
4    0.0190     371.97       138225.9332  86
5    0.0190     361.31       31632.7330   107
6    0.0176     358.13       0.0194       117
7    0.0158     358.11       0.0368       128
8    0.0144     358.12       0.0302       138
9    0.0133     358.09       0.0596       149

--- Optimal solution ---
Number of pipes : 6
Pipe diameter   : 0.0176 m
Max temperature : 358.13 K (84.98 C)
Objective value : 0.0194 K
Feasible        : Yes
Total FEA calls : 149

The optimal configuration here is indicated to have 6 pipes with 17.6 mm diameter, achieving a maximum temperature of 358.13 K, just 0.02 K above the 358.15 K target, well within the permitted 0.1 K tolerance. The full two-sided temperature distribution (exploiting mirroring around the symmetry line for postprocessing) can then be visualized as follows.

postplot( fea, 'surfexpr', 'T' )

fea.grid.p(1,:) = 0.2 - fea.grid.p(1,:);  % Offset and mirror grid points.
postplot( fea_mirror, 'surfexpr', 'T' )

The total optimization process required 149 FEA solves thanks to the combination of decomposition, warm-starting, and caching of FEA results (compared to the several hundred typically needed by a black-box approach).

Running the Example

The complete FEATool simulation script can be downloaded from the link below.

and run directly from the MATLAB command after loading and starting FEATool Multiphysics:

[fea, result] = ex_heattransfer11( 'npipes', [1, 9], 'diam', [0.005, 0.019], 'Tstar', 358.15 );

As shown in this example, the fully programmable simulation API of FEATool Multiphysics makes it straightforward to combine finite element simulations with MATLAB optimization and search routines. The same pattern, wrapping an FEA solve in an objective function and passing it to a minimizer, applies broadly to inverse problems, design optimization, and sensitivity studies across heat transfer, structural mechanics, and fluid dynamics applications. Additional examples of this approach are linked in the sidebar.