Thermal diffusivity: conductivity over heat capacity per volume¶
How fast a temperature spreads through a phase, alpha = k / (rho c_p), from its thermal conductivity and specific heat (lusine.materials.thermal) and its density (mechanical_phases): the coefficient of the heat equation dT/dt = alpha laplacian T.
Source: Incropera & DeWitt, Fundamentals of Heat and Mass Transfer, 6th ed., 2.3 (eq. 2.20). Checked against a reference in numpy, which has its own tests.
- Source
- Incropera & DeWitt, Fundamentals of Heat and Mass Transfer, 6th ed., 2.3 (eq. 2.20).
- Proof
- reference, invariant.
thermal_diffusivity.trait
namespace lusine.properties;
float<m2/s> thermal_diffusivity(float<W/(m K)> conductivity, float<kg/m3> density, float<J/(kg K)> specific_heat)
{
return conductivity / (density * specific_heat);
}
How fast a temperature spreads through a phase, α = k/(ρ c_p), the coefficient of the heat equation ∂T/∂t = α∇²T, from the phase's thermal conductivity and specific heat (lusine.materials.thermal) and its density (mechanical_phases). The library stores no diffusivity: this kernel derives it.
In [1]:
Copied!
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine.library_sources import parse_rows
EPS = float(np.finfo(np.float32).eps) # Trait's float is binary32
The kernel¶
The law as the library states it.
In [2]:
Copied!
print(trait_run.source('properties/thermal_diffusivity.trait'))
print(trait_run.source('properties/thermal_diffusivity.trait'))
// lusine.properties.thermal_diffusivity: how fast a temperature spreads through a phase, alpha = k / (rho c_p), from
// its thermal conductivity and specific heat (lusine.materials.thermal) and its density (mechanical_phases): the
// coefficient of the heat equation dT/dt = alpha laplacian T.
// Source: Incropera & DeWitt, Fundamentals of Heat and Mass Transfer, 6th ed., 2.3 (eq. 2.20).
// Checked against a reference in numpy, which has its own tests.
namespace lusine.properties;
float<m2/s> thermal_diffusivity(float<W/(m K)> conductivity, float<kg/m3> density, float<J/(kg K)> specific_heat)
{
return conductivity / (density * specific_heat);
}
The program¶
Each section runs the kernel over a lattice of items, one a case.
In [3]:
Copied!
print(trait_run.law_program('properties/thermal_diffusivity.trait', 'thermal_diffusivity'))
print(trait_run.law_program('properties/thermal_diffusivity.trait', 'thermal_diffusivity'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<W/(m K)> conductivity;
host float<kg/m3> density;
host float<J/(kg K)> specific_heat;
float<m2/s> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.properties.thermal_diffusivity(items.conductivity, items.density, items.specific_heat); }
Reference: every phase in the library, from its own k, ρ and c_p¶
Each phase with all three, read from the library's tables as committed, is held to k/(ρ c_p) of those numbers, to binary32's rounding. Solid copper's 1.162 × 10⁻⁴ m²/s is within 1% of Incropera & DeWitt's 117 × 10⁻⁶ (Fundamentals of Heat and Mass Transfer, 6th ed., table A.1, pure copper at 300 K, whose ρ is 8933 kg/m³ to the library's 8960); liquid water's 1.455 × 10⁻⁷ m²/s within 1.5% of their 0.147 × 10⁻⁶ (table A.6, at 300 K; the library's water is CoolProp's at 298.15 K, its k 0.6065 W/(m K) to their 0.613).
In [4]:
Copied!
heat = {(r['material'], r['state']): r for r in parse_rows(trait_run.source('materials/generated/thermal.trait'), 'thermal')}
bulk = {(r['material'], r['state']): r for r in parse_rows(trait_run.source('materials/generated/mechanical.trait'), 'mechanical_phases')}
phases = sorted(p for p in heat if p in bulk and 'density' in bulk[p] and 'thermal_conductivity' in heat[p] and 'specific_heat' in heat[p])
k = np.array([heat[p]['thermal_conductivity'] for p in phases])
cp = np.array([heat[p]['specific_heat'] for p in phases])
rho = np.array([bulk[p]['density'] for p in phases])
a = trait_run.law('properties/thermal_diffusivity.trait', 'thermal_diffusivity', conductivity=k, density=rho, specific_heat=cp)
want = k / (rho * cp)
assert np.all(np.abs(a / want - 1) <= 3 * EPS), np.abs(a / want - 1).max()
for p, v in zip(phases, a):
print(f'{p[0]:20s} {p[1]:7s} alpha {v:.4e} m^2/s')
incropera = {('copper', 'solid'): (117e-6, 0.01), ('water', 'liquid'): (0.147e-6, 0.015)}
for p, (v, within) in incropera.items():
got = float(a[phases.index(p)])
assert abs(got / v - 1) < within, (p, got)
print(f"{p[0]} ({p[1]}): {got:.4e} against Incropera & DeWitt's {v:.3e} ({got / v - 1:+.2%})")
heat = {(r['material'], r['state']): r for r in parse_rows(trait_run.source('materials/generated/thermal.trait'), 'thermal')}
bulk = {(r['material'], r['state']): r for r in parse_rows(trait_run.source('materials/generated/mechanical.trait'), 'mechanical_phases')}
phases = sorted(p for p in heat if p in bulk and 'density' in bulk[p] and 'thermal_conductivity' in heat[p] and 'specific_heat' in heat[p])
k = np.array([heat[p]['thermal_conductivity'] for p in phases])
cp = np.array([heat[p]['specific_heat'] for p in phases])
rho = np.array([bulk[p]['density'] for p in phases])
a = trait_run.law('properties/thermal_diffusivity.trait', 'thermal_diffusivity', conductivity=k, density=rho, specific_heat=cp)
want = k / (rho * cp)
assert np.all(np.abs(a / want - 1) <= 3 * EPS), np.abs(a / want - 1).max()
for p, v in zip(phases, a):
print(f'{p[0]:20s} {p[1]:7s} alpha {v:.4e} m^2/s')
incropera = {('copper', 'solid'): (117e-6, 0.01), ('water', 'liquid'): (0.147e-6, 0.015)}
for p, (v, within) in incropera.items():
got = float(a[phases.index(p)])
assert abs(got / v - 1) < within, (p, got)
print(f"{p[0]} ({p[1]}): {got:.4e} against Incropera & DeWitt's {v:.3e} ({got / v - 1:+.2%})")
ash powder alpha 1.3889e-07 m^2/s ceramic liquid alpha 4.6296e-06 m^2/s ceramic solid alpha 3.0166e-07 m^2/s charcoal solid alpha 2.9762e-07 m^2/s clay solid alpha 6.9444e-07 m^2/s coal solid alpha 1.7637e-07 m^2/s coke solid alpha 2.3810e-07 m^2/s copper liquid alpha 4.3549e-05 m^2/s copper solid alpha 1.1625e-04 m^2/s copper_ore solid alpha 1.1023e-06 m^2/s diesel liquid alpha 1.0204e-07 m^2/s dodecane liquid alpha 8.2012e-08 m^2/s graphite solid alpha 3.2050e-06 m^2/s grass solid alpha 1.3889e-07 m^2/s hydrogen liquid alpha 1.4967e-07 m^2/s iron liquid alpha 6.8992e-06 m^2/s iron solid alpha 2.2753e-05 m^2/s iron_ore solid alpha 9.0498e-07 m^2/s kerosene liquid alpha 8.5990e-08 m^2/s leaf solid alpha 3.7037e-07 m^2/s mercury liquid alpha 4.3806e-06 m^2/s mercury solid alpha 4.2047e-06 m^2/s mercury_ore solid alpha 8.5143e-07 m^2/s naphtha liquid alpha 8.4801e-08 m^2/s oil liquid alpha 9.0744e-08 m^2/s rock liquid alpha 3.7135e-07 m^2/s rock solid alpha 8.4388e-07 m^2/s rubber solid alpha 6.5789e-08 m^2/s silicon solid alpha 8.9831e-05 m^2/s silicon_dioxide powder alpha 2.5266e-07 m^2/s soda_lime_glass liquid alpha 5.4348e-07 m^2/s soda_lime_glass solid alpha 4.7619e-07 m^2/s soil solid alpha 4.8077e-07 m^2/s tar solid alpha 6.0606e-08 m^2/s tungsten liquid alpha 1.7623e-05 m^2/s tungsten solid alpha 6.7907e-05 m^2/s tungsten_ore solid alpha 7.6190e-07 m^2/s water liquid alpha 1.4548e-07 m^2/s water solid alpha 1.1966e-06 m^2/s wood solid alpha 1.4706e-07 m^2/s copper (solid): 1.1625e-04 against Incropera & DeWitt's 1.170e-04 (-0.64%) water (liquid): 1.4548e-07 against Incropera & DeWitt's 1.470e-07 (-1.03%)
Invariant: a slab's temperature spreading¶
ρ c_p α is the conductivity again; twice the density halves it and twice the conductivity doubles it, bit for bit. The time a temperature takes to reach a depth L scales as L²/α: across a centimetre of copper it is about a second, across a centimetre of water about eleven minutes.
In [5]:
Copied!
rng = np.random.default_rng(322)
k = np.float32(10 ** rng.uniform(-2, 2.7, 32)).astype(float)
rho = np.float32(10 ** rng.uniform(0, 4.3, 32)).astype(float)
cp = np.float32(10 ** rng.uniform(2, 3.7, 32)).astype(float)
law = lambda **kw: trait_run.law('properties/thermal_diffusivity.trait', 'thermal_diffusivity', **kw)
a = law(conductivity=k, density=rho, specific_heat=cp)
assert np.all(np.abs(a * rho * cp / k - 1) <= 4 * EPS)
assert np.array_equal(law(conductivity=k, density=2 * rho, specific_heat=cp), a / 2)
assert np.array_equal(law(conductivity=2 * k, density=rho, specific_heat=cp), 2 * a)
cu, w = law(conductivity=[401.0, 0.606516], density=[8960.0, 997.048], specific_heat=[385.0, 4181.31])
print(f'L^2 / alpha over 1 cm: copper {1e-4 / cu:.2f} s, water {1e-4 / w / 60:.1f} min')
assert 0.5 < 1e-4 / cu < 1.5 and 10 < 1e-4 / w / 60 < 13
rng = np.random.default_rng(322)
k = np.float32(10 ** rng.uniform(-2, 2.7, 32)).astype(float)
rho = np.float32(10 ** rng.uniform(0, 4.3, 32)).astype(float)
cp = np.float32(10 ** rng.uniform(2, 3.7, 32)).astype(float)
law = lambda **kw: trait_run.law('properties/thermal_diffusivity.trait', 'thermal_diffusivity', **kw)
a = law(conductivity=k, density=rho, specific_heat=cp)
assert np.all(np.abs(a * rho * cp / k - 1) <= 4 * EPS)
assert np.array_equal(law(conductivity=k, density=2 * rho, specific_heat=cp), a / 2)
assert np.array_equal(law(conductivity=2 * k, density=rho, specific_heat=cp), 2 * a)
cu, w = law(conductivity=[401.0, 0.606516], density=[8960.0, 997.048], specific_heat=[385.0, 4181.31])
print(f'L^2 / alpha over 1 cm: copper {1e-4 / cu:.2f} s, water {1e-4 / w / 60:.1f} min')
assert 0.5 < 1e-4 / cu < 1.5 and 10 < 1e-4 / w / 60 < 13
L^2 / alpha over 1 cm: copper 0.86 s, water 11.5 min