Heatsink simulation — Friday, April 30

This version of the program sets the heat-loss parameter ("fraction") to 1/64th, because that value gives a pleasing output. To the extent that it has any physical significance, it suggests that the temperature (and by extension, heat energy) conducts within the heatsink material much more readily than it transfers to the surrounding medium. This is a plausible situation.

heatsink.py

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 30 2021

@author: rmontant
"""
import numpy as np

def pdf(dx):
    return ( np.exp(-(dx**2)/2) / np.sqrt(2*np.pi) )
#--------

class HeatSink:
    surround = -1e100
    initial_temperature = 20

    def __init__(self, nfins=4, finwidth=4, gapwidth=6, baseheight=5):
        self.nfins = nfins
        self.finwidth = finwidth
        self.gapwidth = gapwidth
        self.baseheight = baseheight
        self.ngaps = self.nfins - 1

        ### 2021-04-30 Account for surrounding layers here...
        self.width = nfins * finwidth + self.ngaps * gapwidth + 2
        self.length = self.width
        self.height = int( 1.2 * self.width ) # 1.2 seems like a good number

        ###
        ### 2021-04-22 afternoon:
        ### The cube MUST be initialized to the "surround" value here!
        ### (because it never happens anywhere else...)
        ###
        ### 2021-04-30 ...not here:
        self.cube = np.full(
            (self.width, self.length, self.height), self.surround )

        self.form_base()
        self.form_fins()
        self.build_pointarray()
        self.cube[:,:,0] = 2* self.initial_temperature

        self.tmax = self.initial_temperature
    #--------

    def form_base(self):
        self.cube[1:self.width-1, 1:self.length-1, 1:self.baseheight+1] \
            = self.initial_temperature

    def form_fins(self):
        for f in range(self.nfins):
            fstart = 1 + f * (self.gapwidth + self.finwidth)
            fend = fstart + self.finwidth
            self.cube[ \
                fstart:fend, 1:self.length-1 , self.baseheight+1:self.height-1 \
            ] = self.initial_temperature
    #--------

    def build_pointarray(self):
        coords = []
        for x in range(self.width):
            for y in range(self.length):
                for z in range(self.height):
                    if self.cube[x,y,z] > self.surround:
                        coords.append( (x,y,z) )
        self.coords = np.array( coords )
    #--------

    def set_hotspot(self, w, l, t, scale):
        rmax = ((self.width**2 + self.length**2)**0.5) * scale
        for x in range(1, self.width-1):
            for y in range(1, self.length-1):
                r = ( (w-x)**2 + (l-y)**2 )**0.5
                s = pdf(r/rmax) / pdf(0)

                self.cube[x,y,0] = s * t \
                    + (1-s) * self.cube[x,y,0]
        if t > self.tmax:
            self.tmax = t
    #--------

    def relax_once(self, fraction=(1/8)):
        ###
        ### 2021-04-29:
        ### Make the new cube be a copy of the previous cube - this gets the
        ### hot base set properly, along with the surrounding cells:
        new_cube = self.cube.copy()

        for idx, p in enumerate(self.coords):
            w, l, h = p
            total = []
            for x in (-1,0,1):
                for y in (-1,0,1):
                    for z in (-1,0,1):
                        if self.cube[w+x, l+y, h+z] != self.surround:
                            total.append(self.cube[w+x, l+y, h+z])
                        else:
                            ###
                            ### 2021-04-29:
                            ### Let the edge points leak heat
                            ### into the surrounding points.
                            ### "Fraction" is a heuristic parameter
                            ### that governs the rate of heat loss.
                            ###
                            ### Also: fix a bug:
                            ###     self.cube[w,l,z] --> self.cube[w,l,h]
                            neighbor_temp = \
                                (fraction * self.initial_temperature) \
                                    + ((1 - fraction) * self.cube[w,l,h])
                            total.append(neighbor_temp)
                            #pass

            new_cube[w,l,h] = sum(total)/len(total)
        return new_cube
    #--------

#--------

model-heatsink.py

Two fins, "fraction" settable from the command line or user input.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
2021-04-20 "final" version
@author: rmontant
"""
import sys, time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize
from mpl_toolkits.mplot3d import Axes3D

import heatsink

def draw_contour(axis, surface, tmin=None, tmax=None):
    if not tmin:
        tmin = np.min(surface)
    if not tmax:
        tmax = np.max(surface)
    normal = Normalize(tmin, tmax)
    pcolormap = axis.pcolormesh(surface, cmap=cm.afmhot, norm=normal)
    axis.contour(surface, cmap=cm.hot)
    plt.colorbar(pcolormap, ax=axis)
#--------

def rgb(v):
    return (np.sqrt(v), v, 0.9*v, 0.05)
#--------

def draw_heatsink(axis, heatsink):
    xcoords = [ p[0] for p in heatsink.coords]
    ycoords = [ p[1] for p in heatsink.coords]
    zcoords = [ p[2] for p in heatsink.coords]
    colors = [ rgb(0.6 * z / heatsink.height) for z in zcoords]
    axis.scatter(xcoords, ycoords, zcoords, color=colors, marker=',')
    axis.set_xlabel('width')
    axis.set_ylabel('length')
    axis.set_zlabel('height')
#--------

###
### 2021-04-30 --- add a name parameter, to distinguish the saved plots.
###
def display_results(heatsink, name):
    figA = plt.figure(figsize=(10,10))
    figA.suptitle(name+'-z')
    for layer in range(9):
        layer_idx = round((layer / 9) * heatsink.height)
        surface = heatsink.cube[ : , : , layer_idx]
        axis = figA.add_subplot(3,3, (layer+1))
        draw_contour(axis, surface,
                    ### Force the temperature scale to something impressive:
                     tmin=heatsink.initial_temperature,
                     #tmin=100,
                     tmax=heatsink.tmax)
    plt.savefig(name+'-z.png')

    ###
    ### 2021-4-30 --- add slicing along the length,
    ### to show the fins in cross-section:
    figB = plt.figure(figsize=(10,10))
    figB.suptitle(name+'-l')
    for layer in range(9):
        layer_idx = round((layer / 9) * heatsink.length)
        surface = heatsink.cube[ : , layer_idx, : ]
        axis = figB.add_subplot(3,3, (layer+1))
        draw_contour(axis, surface,
                    ### Force the temperature scale to something impressive:
                     tmin=heatsink.initial_temperature,
                     #tmin=100,
                     tmax=heatsink.tmax)
    plt.savefig(name+'-l.png')
    plt.show()
#--------


def main(argv=[__name__]):
    if '-f' in argv:
        idx = argv.index('-f')
        loss_fraction = int(argv[idx+1])
    else:
        loss_fraction = int(input('heat-loss fraction? '))

    htsnk = heatsink.HeatSink(nfins=2, finwidth=3, gapwidth=4, baseheight=4)

    id_string = 'heatsink-r.{:d}-{:d}-{:d}-{:d}.{:d}' \
        .format(htsnk.nfins, htsnk.finwidth, htsnk.gapwidth, htsnk.baseheight, loss_fraction)

    nhotspots = 3
    for i in range(nhotspots):
        temperature = 100 + 10 * np.random.random()
        w = np.random.randint(1, htsnk.width-1)
        l = np.random.randint(1, htsnk.length-1)
        scale = 0.07 + 0.06 * np.random.random()
        print('hotspot {:d} is {:.3} at ({:d}, {:d})' \
              .format(i, scale, w, l))
        htsnk.set_hotspot(w, l, temperature, scale)

    #----
    fig = plt.figure(figsize=(8,8))
    fig.suptitle(id_string)
    a1 = fig.add_subplot(2, 2, 1)
    draw_contour(a1, htsnk.cube[:,:,0], \
                 tmin=htsnk.initial_temperature,
                 tmax=htsnk.tmax
    )

    a2 = fig.add_subplot(2,2, 2)
    surface = htsnk.cube[:, htsnk.length//2, :]
    draw_contour(a2, surface,
                 tmin=htsnk.initial_temperature,
                 tmax=htsnk.tmax,
    )

    a4 = fig.add_subplot(2,2, 4, projection='3d')
    draw_heatsink(a4, htsnk)

    plt.savefig(id_string + '.png')
    #----

    display_results(htsnk, id_string)


    timestart = time.process_time()
    changed = True
    passes = 0
    while changed:
        passes += 1
        elapsed = time.process_time() - timestart
        min = int(elapsed / 60)
        sec = elapsed % 60
        print('\r{:5d} passes, {:d}m {:.3f}s'.format(passes, min,sec),end='')

        ### 2021-04-30 --- slower heat absorption
        new_cube = np.round( htsnk.relax_once(fraction=(1/loss_fraction)), 3 )
        changed = np.any( new_cube != htsnk.cube )
        # or...
        # changed = not np.all( new_heatsink.cube == htsnk.cube )
        htsnk.cube = new_cube


    ###
    ### 2021-04-30 --- add timestamps to output
    ###
    timedone = time.process_time()
    elapsed = timedone - timestart
    min = int(elapsed / 60)
    sec = elapsed % 60
    print('\r{:5d} passes total, {:d}m {:.3f}s elapsed'.format(passes, min,sec))

    display_results(htsnk, id_string)
#--------

if __name__ == '__main__':
    sys.exit(main(sys.argv))
#--------

Output — "fraction" = 1/4

$> python3 model-heatsink.py
hotspot 0 is 0.0786 at (10, 2)
hotspot 1 is 0.11 at (8, 6)
hotspot 2 is 0.0708 at (3, 10)
 8115 passes total, 53m 25.838s elapsed
$> 
initial view   cross-sections along z-axis   cross-sections along y-axis