## 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.

Geometry

• 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

• 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);
```

Create table

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

Import geometries into the PDE models

```for k = 1:size(T,1)
importGeometry(T.Model{k},T.Properties.RowNames{k});
% 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');
end
end
```

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

```figure
pdegplot(T.Model{1},'FaceLabels','on');
title(T.Properties.RowNames{1});
view(0,90);
figure
pdegplot(T.Model{end},'FaceLabels','on');
title(T.Properties.RowNames{end});
view(0,90);
```  ## 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);
```

```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
model.generateMesh('hmax',0.25/4);
[p,e,t] = meshToPet(model.Mesh);
% variable component of boundary conditions is generalized Neumann BC
% on maxMinSpreadFaceNodes and also face on center protrusion
'q',hc);
% 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 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
res(find(isnan(res)))=0;
return
end
area = 2*offsetX*(maxZ-minZ);
avgTemp = integral2(@intFun,-offsetX,+offsetX,minZ,maxZ)/area;
end
```

## 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));
end
```

Populate table with results of computations that are performed asynchronously

```for idx = 1:size(T,1)
% get result for next geometry that was solved
% 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);
end
```
```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');
TOpCost(1:5,:)
```
```ans =

__________________    ________    ___________    _______________    __________________    ____________    _____________

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;
figure
pdeplot3D(TOpCost.Model{1},'colormapdata',u_Optimal);
view(45,90);
snapnow
view(114,51);
snapnow
```  Top-5 max-min spread sorted from smallest to largest

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

__________________    ________    ___________    _______________    __________________    ____________    _____________

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));
figure
view(45,90);
snapnow
view(114,51);
snapnow
```  ## Takeaways

• 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
```end
```