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.
Step.1: Adding thermal conductivities
First, the elastic material properties (lines 8-10) are changed to the thermal conductivities of materials
k0 = 1; % Good thermal conductivity kmin = 1e-3; % Poor thermal conductivity
Step.2: Changing Boundary conditions
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).
% 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);
Step.3: Changing auxiliary variables
Since there is only one DOF per node in heat condition problems, some variables need to change correspondingly, such as , .
Change the total number of DOFs and the force vector in lines 21-22
ndof = (nelx+1)*(nely+1)*(nelz+1); F = sparse(1:ndof,1,-0.01,ndof,1);
The element conductivity matrix is called in line 25 by
KE = lk_H8(k0);
and it is defined in lines 99-145
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
The finite element connectivity matrix is changed in lines 30-35
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*8*nele,1); jK = reshape(kron(edofMat,ones(1,8))',8*8*nele,1);
Step.4: Changing Finite Element Analysis Functions
The global conductivity matrix is assembled in a different way, hence line 70 need to change as
sK = reshape(KE(:)*(kmin+(1-kmin)*xPhys(:)'.^penal),8*8*nele,1);
The evulation of the objective function and analysis of the sensitivity are given in lines 75-76
c = sum(sum(sum((kmin+(1-kmin)*xPhys.^penal).*ce))); dc = -penal*(1-kmin)*xPhys.^(penal-1).*ce;
Step.5: Running program
The optimized topology is derived as shown below by promoting the following line in the :