MATLAB Examples


Best Placement of Holes in Cylinder to Achieve Target Average Temperature - A Parametric Study

This examples conducts a parametric study in which heat conduction simulation is performed over a set of similar geometries to determine which geometry "best" meets an average temperature on an specified output area. The geometry is a cylinder like structure as shown below and has a ring of holes running longitudinally through the structure. The problem has the following characteristics:

Boundary conditions

  • The input heat source is applied on the faces of the holes.
  • The longitudinal surface and the surface on the center protrusion have convective boundary conditions. The output surface is the rectangular subsurface of the longitudinal surface. The objective of the example is to achieve a target average temperature on output surface in the "best" manner possible as described later.
  • All other faces not indicated above are insulated and thus have zero Neumann boundary conditions.


  • Each geometry has a unique pair of (#holes, radius of ring of holes). All other geometry parameters are held constant.
  • Although it is possible to exploit symmetry in some of the geometries in order to reduce the problem to 2 geometry dimensions, it was not done so in this example.


  • Results are collected from all simulations and the best geometry in terms of lowest max-min temperature spread on the longitudinal face and the best geometry for lowest operating cost (input flux) are identified.
  • The implementation uses 'parfeval' from the Parallel Computing Toolbox to speed up the parametric study.
function heating

Import candidate geometries and create models

The STL files are read - each file corresponds to a different parameter-pair: (#holes, radius of ring of holes).

fileList = ls('cyl_*.STL');
fileList = mat2cell(fileList,ones(size(fileList,1),1));

The PDE is a scalar, laplace equation

N = 1;

A table will be created to organize data and results from all runs

Table column corresponding to the PDE models

Model = cellfun(@(~) createpde(N),fileList,'UniformOutput',false);

Table columns corresponding to #holes and radius (note: this is the radius of the ring and not the radius of the hole) are extracted.

paramList = cellfun(@(fileName) regexpi(fileName,'cyl_(.*)_(.*).STL','tokens'),fileList, 'UniformOutput',false);
NumHoles = cellfun(@(entry) str2double(entry{1}(1)),paramList,'UniformOutput',false);
HolesRadius = cellfun(@(entry) str2double(entry{1}(2)),paramList,'UniformOutput',false);

Create table

T = [table(Model), cell2table(NumHoles,'RowNames',fileList), cell2table(HolesRadius)];

Import geometries into the PDE models

for k = 1:size(T,1)
    % The relation of faces to holes is known; report errors for unexpected
    % relation
    if T.Model{k}.Geometry.NumFaces ~= (3 + T.NumHoles(k) + 2)
        error('unexpected number of faces');

Sort table, first by #holes and then by radius

T = sortrows(T,{'NumHoles','HolesRadius'},{'ascend','ascend'});

Plot two extreme geometries to show the range of geometry variations


Input setup

Ambient temperature

ambientTemp = 6;

Target average nodal temperature on output surface

targetTemp = 15;

PDE coefficients for laplace equation (heat conduction)

c = 1e-1;
a = 0;
f = 0;

Any face in (inputFacesBegin:(inputFacesBegin +numHoles)) is an input heat source face

inputFacesBegin = 4;

Output setup

Table column for capturing max-min temperature on output area; it is desirable to have a low spread

MaxMinSpread = zeros(size(T,1),1);

Table column for operating cost (total flux going into solid via the input heat source faces); it is desirable to minimize this

OperatingCost = zeros(size(T,1),1);

AvgTempVariable column corresponds to the variable contribution towards the average temperature solution on the output surface. The constant contribution is simply ambientTemp.

AvgTempVariable = zeros(size(T,1),1);

Table column for scale factor for AvgTempVariable to help match targetTemp

InputForTargetTemp = zeros(size(T,1),1);

Add these columns to table

T = [T table(AvgTempVariable,InputForTargetTemp,MaxMinSpread,OperatingCost)];

Face on which the max-min temperature spread is measured.

MaxMinSpreadFace = 1;

Output surface in XZ plane: -offsetX:offsetX, offsetY, minZ:maxZ where average temperature is calculated

offsetY = -1.875;
offsetX = sqrt(2^2-offsetY^2);
minZ = 0;
maxZ = 1;

Convective heat transfer coefficient

hc = 0.3;

Function for applying boundary conditions, meshing, and solving per geometry

    function [avgTemp,maxMinSpread,resultVariableBC] = solveGeometry(model,numHoles)
        % generate mesh with 'hmax' = 1/4th of hole radius
        % extract nodes on MaxMinSpreadFace
        [p,e,t] = meshToPet(model.Mesh);
        maxMinSpreadFaceNodes = e.getNodes(MaxMinSpreadFace);
        % variable component of boundary conditions is generalized Neumann BC
        % on maxMinSpreadFaceNodes and also face on center protrusion
        % apply unit flux on input heat source faces
        model.applyBoundaryCondition('Face',(inputFacesBegin:(inputFacesBegin + numHoles)),'g',1);
        % solve to get result for variable BC
        resultVariableBC = assempde(model,c,a,f);
        % calculate max-min spread
        t1 = resultVariableBC(maxMinSpreadFaceNodes);
        maxMinSpread = max(t1) - min(t1);
        % calculate average temp. on rectangular output surface area
        myInterpolant = pdeInterpolant(p,t,resultVariableBC);
        function res = intFun(x,z)
            % Y offset is pushed a bit inwards to avoid missing data along
            % Y axis
            res = evaluate(myInterpolant,x,(offsetY+0.01)*ones(size(x)),z);
            res = reshape(res,size(x));
            % NaNs dues to XZ plane overshoot are set to zero
        area = 2*offsetX*(maxZ-minZ);
        avgTemp = integral2(@intFun,-offsetX,+offsetX,minZ,maxZ)/area;

Solve for all geometries

Use parfeval to perform asynchronous computation and speed up overall simulation time

pool = gcp();
Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers.

Futures are created for the geometries

for idx = 1:size(T,1)
    % coefficients of PDE
    F(idx) = parfeval(pool,@solveGeometry,3,T.Model{idx},T.NumHoles(idx));

Populate table with results of computations that are performed asynchronously

for idx = 1:size(T,1)
    % get result for next geometry that was solved
    [completedIdx,avgTemp,maxMinSpread] = fetchNext(F);
    % Set the average temperature contribution of the variable part of the
    % BC
    T.AvgTempVariable(completedIdx) = avgTemp;
    fprintf('%d of %d models simulated\n',completedIdx,size(T,1));
    % compute scale factor for the contribution of the variable part. As
    % mentioned earlier ambientTemp corresponds to contribution of the constant part.
    T.InputForTargetTemp(completedIdx) = (targetTemp-ambientTemp)./T.AvgTempVariable(completedIdx);
    T.MaxMinSpread(completedIdx) = maxMinSpread;
2 of 41 models simulated
1 of 41 models simulated
3 of 41 models simulated
4 of 41 models simulated
5 of 41 models simulated
6 of 41 models simulated
7 of 41 models simulated
8 of 41 models simulated
9 of 41 models simulated
10 of 41 models simulated
11 of 41 models simulated
12 of 41 models simulated
13 of 41 models simulated
14 of 41 models simulated
15 of 41 models simulated
16 of 41 models simulated
17 of 41 models simulated
18 of 41 models simulated
20 of 41 models simulated
19 of 41 models simulated
21 of 41 models simulated
22 of 41 models simulated
23 of 41 models simulated
24 of 41 models simulated
25 of 41 models simulated
26 of 41 models simulated
27 of 41 models simulated
28 of 41 models simulated
29 of 41 models simulated
30 of 41 models simulated
31 of 41 models simulated
32 of 41 models simulated
33 of 41 models simulated
34 of 41 models simulated
36 of 41 models simulated
37 of 41 models simulated
35 of 41 models simulated
39 of 41 models simulated
40 of 41 models simulated
38 of 41 models simulated
41 of 41 models simulated

Report and visualize results

Calculate operating cost

T.OperatingCost = T.InputForTargetTemp.*T.NumHoles;

Top-5 operating cost sorted from smallest to largest

TOpCost = sortrows(T,'OperatingCost');
ans = 

                           Model           NumHoles    HolesRadius    AvgTempVariable    InputForTargetTemp    MaxMinSpread    OperatingCost
                     __________________    ________    ___________    _______________    __________________    ____________    _____________

    cyl_4_1.2.STL    [1x1 pde.PDEModel]    4           1.2            1.9545             4.6048                 0.7581         18.419       
    cyl_4_1.1.STL    [1x1 pde.PDEModel]    4           1.1            1.8816             4.7831                0.60484         19.133       
    cyl_4_1.0.STL    [1x1 pde.PDEModel]    4             1            1.8197             4.9458                0.49596         19.783       
    cyl_3_1.2.STL    [1x1 pde.PDEModel]    3           1.2            1.3607              6.614                0.69789         19.842       
    cyl_3_1.1.STL    [1x1 pde.PDEModel]    3           1.1             1.354             6.6468                0.56056          19.94       

Plot result for geometry with lowest operating cost

[~,~,resultVariableBC] = solveGeometry(TOpCost.Model{1},TOpCost.NumHoles(1));
u_Optimal = resultVariableBC*TOpCost.InputForTargetTemp(1) + ambientTemp;

Top-5 max-min spread sorted from smallest to largest

TMaxMinSpread = sortrows(T,'MaxMinSpread');
ans = 

                           Model           NumHoles    HolesRadius    AvgTempVariable    InputForTargetTemp    MaxMinSpread    OperatingCost
                     __________________    ________    ___________    _______________    __________________    ____________    _____________

    cyl_3_0.9.STL    [1x1 pde.PDEModel]    3           0.9             1.3254            6.7903                0.41703         20.371       
    cyl_4_0.9.STL    [1x1 pde.PDEModel]    4           0.9             1.7654             5.098                0.44566         20.392       
    cyl_3_1.0.STL    [1x1 pde.PDEModel]    3             1              1.342            6.7062                 0.4798         20.119       
    cyl_2_0.9.STL    [1x1 pde.PDEModel]    2           0.9            0.71197            12.641                0.48223         25.282       
    cyl_5_0.9.STL    [1x1 pde.PDEModel]    5           0.9             2.1201             4.245                0.48678         21.225       

Plot result for geometry with lowest max-min spread

[~,~,resultVariableBC] = solveGeometry(TMaxMinSpread.Model{1},TMaxMinSpread.NumHoles(1));
u_Optimal = resultVariableBC*TMaxMinSpread.InputForTargetTemp(1) + ambientTemp;


  • Programmatically solve PDEs and perform parametric studies. Useful if there is little fundamental variation between different design points.
  • Perform custom post-processing with useful reports
  • Use Parallel Computing Toolbox to accelerate simulations