Hilbert curve heater track
A heater track along a Hilbert curve fills its square evenly and therefore heats it evenly. The example has two parts. The first builds the track of order 7, with 16384 points on a pitch = 2 mm grid and width = 1 mm wide. The second puts a heater of order 4 on an aluminum plate and computes the current through the track, its Joule heat and the plate’s temperature.
[1]:
import math
import numpy as np
import ngsolve as ngs
from ngsolve_webgpu.jupyter import Draw
import planex as px
hilbert(order) gives the points of the curve, on the grid of the pitch:
[3]:
pitch = 2.0
width = 1.0
def hilbert(order):
n, points = 2**order, []
for d in range(n * n):
x = y = 0
s, t = 1, d
while s < n:
rx = 1 & (t // 2)
ry = 1 & (t ^ rx)
if ry == 0:
if rx == 1:
x, y = s - 1 - x, s - 1 - y
x, y = y, x
x, y = x + s * rx, y + s * ry
t //= 4
s *= 2
points.append((pitch * x, pitch * y))
return points
points = hilbert(7)
print(len(points), "points")
16384 points
The track
The path is a sketch of lines through the points, and the track is its wire offset by half the width:
[4]:
path = px.QuickSketch(points[0])
for p in points[1:]:
path.Line(p)
track = path.Wire().Offset(width / 2)
Check: the track’s area is the width times the path length, plus the two round ends. At every 90° turn, the square (w/2)² inside the corner is counted twice, and the quarter disc outside it is added.
[5]:
turns = sum(1 for p, q, r in zip(points, points[1:], points[2:])
if (q[0] - p[0]) * (r[1] - q[1]) != (q[1] - p[1]) * (r[0] - q[0]))
length = pitch * (len(points) - 1)
exact = width * length + math.pi * width**2 / 4 + turns * (width / 2) ** 2 * (math.pi / 4 - 1)
print(turns, "turns")
print("area", track.mass, "exact", exact)
13106 turns
area 32063.642480531933 exact 32063.642480535134
A heater on a plate
The track of order 7 is 33 m long, far too long to mesh here. The heater is the same curve of order 4 on the same grid, with 256 points in a 30 mm square, made of constantan foil 25 µm thick, on a 40 mm aluminum plate 1.5 mm thick. Leads 3 long take both ends to the plate’s lower edge, 2 below the heater, whose two halves are named as the contacts. The plate reaches 5 beyond the heater on the other sides. Cutting the track to the plate gives its ends those names, and gluing it into the plate gives one mesh for both (Figure 1).
[6]:
heater_points = hilbert(4)
side = heater_points[-1][0]
path = px.QuickSketch((0, -3)).Line(heater_points[0])
for p in heater_points[1:]:
path.Line(p)
path.Line((side, -3))
plate = px.QuickSketch((-5, -2)).Line((side / 2, -2), name="minus").Line((side + 5, -2), name="plus")
plate = plate.Line((side + 5, side + 5)).Line((-5, side + 5)).Close().Face()
plate.name = "plate"
heater = path.Wire().Offset(width / 2) * plate
heater.name = "track"
heater.vertices.hpref = 1
board = px.Glue([plate, heater])
Current
In the foil the potential solves Laplace’s equation with the supply voltage on one contact and zero on the other. The current comes from the leftover equations at the fixed contact, so the electric power V·I must equal the Joule heat integrated over the track. The resistance in squares (the resistance over the foil’s sheet resistance ρ/t) is less than the center line’s length over the width, because the current cuts across at a corner. A sharp 90° corner counts as about 0.56 of a square instead of one, and these corners have a rounded outside. The inside corners of the turns are re-entrant, so the current density is singular there, and the mesh is refined towards them, with the potential of order
The color scale ends at 1.2 times the mean field along the track, and at the inside corners the field goes far above it.
[8]:
rho = 4.9e-4
foil = 0.025
volts = 5.0
mesh = ngs.Mesh(px.OCCGeometry(board).GenerateMesh(maxh=0.5))
mesh.RefineHP(levels=2)
mesh.Curve(3)
E = ngs.H1(mesh, order=3, definedon="track", dirichlet="plus|minus")
phi, psi = E.TnT()
conduct = ngs.BilinearForm(foil / rho * ngs.grad(phi) * ngs.grad(psi) * ngs.dx("track", bonus_intorder=4)).Assemble()
potential = ngs.GridFunction(E)
potential.Set(volts, definedon=mesh.Boundaries("plus"))
potential.vec.data += conduct.mat.Inverse(E.FreeDofs()) * (-conduct.mat * potential.vec)
reaction = potential.vec.CreateVector()
reaction.data = conduct.mat * potential.vec
contact = ngs.GridFunction(E)
contact.Set(1, definedon=mesh.Boundaries("plus"))
current = ngs.InnerProduct(contact.vec, reaction)
joule = foil / rho * ngs.grad(potential) * ngs.grad(potential)
center_line = [(0, -2)] + heater_points + [(side, -2)]
corners = sum(1 for p, q, r in zip(center_line, center_line[1:], center_line[2:])
if (q[0] - p[0]) * (r[1] - q[1]) != (q[1] - p[1]) * (r[0] - q[0]))
strip = sum(math.dist(p, q) for p, q in zip(center_line, center_line[1:])) / width
squares = volts / current * foil / rho
heat = ngs.Integrate(joule, mesh, definedon=mesh.Materials("track"), order=10)
print(mesh.ne, "elements")
print("current", round(current, 4), "A, resistance", round(volts / current, 3), "Ohm")
print("power V I", volts * current, "W, Joule heat", heat, "W")
print(round(squares, 1), "squares instead of", round(strip), "along the center line")
print("each of the", corners, "corners counts", round(1 - (strip - squares) / corners, 2), "of a square")
Draw(-ngs.grad(potential), mesh, min=0, max=1.2 * volts / (squares * width), autoscale=False)
24080 elements
current 0.6025 A, resistance 8.298 Ohm
power V I 3.0126088917675964 W, Joule heat 3.012608891758272 W
423.4 squares instead of 514 along the center line
each of the 206 corners counts 0.56 of a square
Temperature
The plate conducts the heat sideways and gives it to still air on both faces. With the air at 20 °C, the heat the plate gives off must equal the electric power; the refined elements at the corners are not affine, so the heat source is integrated with the same rule as the Joule heat above. The temperature over the heated square shows how even the heater is.
[9]:
k = 0.237
plate_t = 1.5
h = 1e-5
T_air = 20.0
T = ngs.H1(mesh, order=2)
theta, v = T.TnT()
plate_form = k * plate_t * ngs.grad(theta) * ngs.grad(v) + 2 * h * theta * v
heat = ngs.BilinearForm(plate_form * ngs.dx(bonus_intorder=4)).Assemble()
rules = {ngs.ET.TRIG: ngs.IntegrationRule(ngs.TRIG, 10), ngs.ET.QUAD: ngs.IntegrationRule(ngs.QUAD, 10)}
source = ngs.LinearForm(joule * v * ngs.dx("track", intrules=rules)).Assemble()
rise = ngs.GridFunction(T)
rise.vec.data = heat.mat.Inverse(T.FreeDofs()) * source.vec
print("heat given off", ngs.Integrate(2 * h * rise, mesh, order=10), "W, electric power", volts * current, "W")
points3 = {ngs.ET.TRIG: ngs.IntegrationRule(ngs.TRIG, 3), ngs.ET.QUAD: ngs.IntegrationRule(ngs.QUAD, 3)}
samples = mesh.MapToAllElements(points3, ngs.VOL)
xy = ngs.CF((ngs.x, ngs.y))(samples)
inside = (xy[:, 0] > 0) & (xy[:, 0] < side) & (xy[:, 1] > 0) & (xy[:, 1] < side)
temperature = T_air + np.asarray(rise(samples)).ravel()[inside]
print("heated square from", round(temperature.min(), 2), "to", round(temperature.max(), 2), "°C")
Draw(T_air + rise, mesh)
heat given off 3.0126088917606317 W, electric power 3.0126088917675964 W
heated square from 121.56 to 121.98 °C
