Heat in physics is defined as energy transferred between a system and its surrounding. The direct microscopic exchange of kinetic energy of particles through the boundary between two systems is called diffusion or heat conduction. In this tutorial, you will learn how to use top3d to solve heat conduction problems.
The implementation of heat conduction problems is not more complex than the one for compliant mechanism synthesis since the number of Degree Of Freedom (DOF) per node is one rather than three. Following the implementation of heat conduction problems in two dimensions, the implementation for three dimension problem as shown below is suggested in the following steps.

[expand title="Step.1: Adding thermal conductivities" expanded="true" tag="h5"]

First, the elastic material properties (lines 8-10) are changed to the thermal conductivities of materials

[matlab firstline="8"]
k0 = 1; % Good thermal conductivity
kmin = 1e-3; % Poor thermal conductivity
[/matlab]

[/expand]

[expand title="Step.2: Changing Boundary conditions" expanded="true" tag="h5"]

The boundary conditions for the heat condition problem, i.e., a rectangular plate with a heat sink on the middle of top face and all nodes are given a thermal load as shown above, are changed corresponding (lines 10-18).

[matlab firstline="10"]
% USER-DEFINED SUPPORT FIXED DOFs
il = nelx/2-nelx/20:nelx/2+nelx/20; jl = nely; kl = 0:nelz;
fixedxy = il*(nely+1)+(nely+1-jl);
fixednid = repmat(fixedxy',size(kl)) + ...
repmat(kl*(nelx+1)*(nely+1),size(fixedxy,2),1);
fixeddof = reshape(fixednid,[],1);
[/matlab]

[/expand]

[expand title="Step.3: Changing auxiliary variables" expanded="true" tag="h5"]

Since there is only one DOF per node in heat condition problems, some variables need to change correspondingly, such as \$latex \texttt{ndof}\$, \$latex \texttt{edofMat}\$.

Change the total number of DOFs and the force vector in lines 21-22

[matlab firstline="21"]
ndof = (nelx+1)(nely+1)(nelz+1);
F = sparse(1:ndof,1,-0.01,ndof,1);
[/matlab]

The element conductivity matrix is called in line 25 by

[matlab firstline="25"]KE = lk_H8(k0);[/matlab]

and it is defined in lines 99-145

[matlab firstline="99"]
function [KE] = lk_H8(k)
A1 = 4*eye(2); A2 = -eye(2);
A3 = fliplr(A2); A4 = -ones(2);
KE1 = [A1 A2; A2 A1];
KE2 = [A3 A4; A4 A3];
KE = 1/12 * k * [KE1 KE2; KE2 KE1];
end
[/matlab]

The finite element connectivity matrix \$latex \texttt{edofMat}\$ is changed in lines 30-35

[matlab firstline="30"]
edofVec = nodeids(:)+1;
edofMat = repmat(edofVec,1,8)+ ...
repmat([0 nely + [1 0] -1 ...
(nely+1)(nelx+1)+[0 nely + [1 0] -1]],nele,1);
iK = reshape(kron(edofMat,ones(8,1))',8
8nele,1);
jK = reshape(kron(edofMat,ones(1,8))',8
8*nele,1);
[/matlab]

[/expand]

[expand title="Step.4: Changing Finite Element Analysis Functions" expanded="true" tag="h5"]

The global conductivity matrix is assembled in a different way, hence line 70 need to change as

[matlab firstline="70"]sK = reshape(KE(:)(kmin+(1-kmin)xPhys(:)'.^penal),88nele,1);[/matlab]

The evulation of the objective function and analysis of the sensitivity are given in lines 75-76

[matlab firstline="75"]
c = sum(sum(sum((kmin+(1-kmin)*xPhys.^penal).ce)));
dc = -penal
(1-kmin)*xPhys.^(penal-1).*ce;
[/matlab]

[/expand]

[expand title="Step.5: Running program" expanded="true" tag="h5"]

The optimized topology is derived as shown below by promoting the following line in the \$latex \textsc{Matlab}\$:

[matlab light="true"]top3d(40,40,5,0.30,3.0,1.4)[/matlab]

[/expand]