Skip to content

Repository files navigation

CryoHeatFlow

A Python package for cryogenic thermal analysis and heat transfer calculations. This package provides functions for calculating thermal conductivity, thermal power transfer, thermal boundary conductance, and multilayer insulation effectiveness.

image

Installation

pip install cryoheatflow

Table of Contents

Features

  • Thermal conductivity calculations for various materials at cryogenic temperatures
  • Thermal power transfer through conductors and insulators
  • Thermal boundary conductance across joints and interfaces
  • Multilayer insulation effectiveness calculations
  • Coax cable heat flow simulation with the outer conductor anchored at multiple temperature stages
  • Area calculations for various cross-sectional geometries (coax, AWG wire, etc)

Quick Start

Calculate Thermal Conductivity

import cryoheatflow

# Get thermal conductivity for stainless steel at 10K
k_conductivity_function = cryoheatflow.conductivity.k_ss
T = 10  # Temperature in Kelvin
result = k_conductivity_function(T)
print(f'Thermal conductivity = {result} W/m*K')

Calculate Thermal Power Transfer

Let's say you wanted to connect a stainless-steel microwave coax line from a 40K stage to a 4K stage. The coax has a diameter of 0.085" (so-called "085" coax), and is 30mm long. How much heat would be transferred?

import cryoheatflow

# Select stainless steel as the material
k = cryoheatflow.conductivity.k_ss
area = cryoheatflow.area.coax_085  # 0.085" outer-diameter coax
length = 30e-3  # 30 mm
T1 = 40  # 40 K 
T2 = 4   # 4 K

P, G, R = cryoheatflow.calculate_thermal_transfer(k, area, length, T1, T2)
print(f'Power transmission = {P*1e3:0.3f} mW')
print(f'Thermal conductance = {G:0.6f} W/K')
print(f'Thermal resistance = {R:0.3f} K/W')

giving us

Power transmission = 4.844 mW
Thermal conductance = 0.000135 W/K
Thermal resistance = 7432.015 K/W

Calculate thermal boundary conductance

Now let's say you want to anchor a 1.5x1.5 cm^2 sample to your 4K stage, and you put grease between the sample and the 4K stage. Your sample is going to generate 2 mW of heat load and going to warm up a little. What temperature is your sample going to be at?

First, we calculate the thermal boundary conductance (in watts per kelvin), and/or its inverse quantity, the thermal boundary resistance:

import cryoheatflow

T = 4  # Temperature in Kelvin
area_m2 = 15e-3 * 15e-3  # 15 mm x 15 mm contact area

h = cryoheatflow.conductivity.h_grease(T=T, area=area_m2)
print(f'Thermal conductance = {h:0.3f} W/K')
print(f'Thermal resistance = {1/h:0.3f} K/W')

This gives us Thermal resistance R = 38.384 K/W. We can then estimate the temperature by the simple relation

(temperature increase) = (thermal resistance) x (heating power)

Giving us a temperature increase of ~76.8 mK.

Calculate Temperature Rise

If you have a thermal conductor with a known heat load applied at one end and the other end anchored at a known temperature, you can calculate the temperature rise at the hot end.

For example, assume you have a 4mm thick x 2mm wide x 100mm long strip of aluminum 6061-T6 that's attached to a 40K coldhead at one end. If the other end of the strip has 250 mW of heat load applied to it, what will the temperature be at the hot end?

import cryoheatflow

k = cryoheatflow.conductivity.k_al6061
area = 4e-3 * 2e-3  # 4 mm x 2 mm
length = 100e-3  # 100 mm
T1 = 40  # 40 K (cold end temperature)
heat_load = 0.25  # 250 mW

T2, thermal_conductance, thermal_resistance = cryoheatflow.calculate_temperature_rise(k, area, length, T1, heat_load)
print(f'Temperature at cold end = {T1:0.3f} K')
print(f'Temperature at hot end = {T2:0.3f} K')
print(f'Thermal conductance = {thermal_conductance:0.3f} W/K')
print(f'Thermal resistance = {thermal_resistance:0.3f} K/W')

giving us

Temperature at cold end = 40.000 K
Temperature at hot end = 83.680 K
Thermal conductance = 0.006 W/K
Thermal resistance = 174.719 K/W

Multilayer Insulation Analysis

from cryoheatflow import solve_multilayer_insulation
from cryoheatflow.emissivity import Al_polished, Al_oxidized, mylar

# Calculate effectiveness of multilayer insulation
T1 = 4   # Cold side temperature (K)
T2 = 85  # Warm side temperature (K)
N = 2    # Number of mylar layers
emissivity1 = Al_oxidized    # Emissivity of the first layer (e.g. 300K walls)
emissivity_mylar = mylar     # Emissivity of the multilayer mylar layers
emissivity2 = Al_polished    # Emissivity of the last layer (e.g. 40K walls)
area = (20e-2)**2           # Area in m^2

layer_temps, qdot = solve_multilayer_insulation(T1, T2, N, emissivity1, emissivity_mylar, emissivity2, area)
print(f'Layer temperatures: {layer_temps}')
print(f'Thermal power: {abs(qdot)} W')

Coax Cable Heat Flow Simulation

When running a coax cable from room temperature down through a cryostat, you can typically only heat-sink the outer conductor at the intermediate stages — the inner conductor is only thermalized at the two ends (e.g. at the 300K feedthrough and at the 4K device). Heat then sneaks down the inner conductor, partially thermalizing to the outer conductor through the dielectric along the way. solve_coupled_conductor_heat_flow() solves for the temperature profiles of both conductors and reports how much heat is deposited into each temperature stage.

import numpy as np
from cryoheatflow import solve_coupled_conductor_heat_flow, k_ss, k_ptfe

# 085 stainless coax with PTFE dielectric, 10 cm long
in2m = 1/39.37
d_inner = 0.020*in2m       # inner conductor diameter
d_dielectric = 0.066*in2m  # dielectric outer diameter
d_outer = 0.085*in2m       # outer conductor outer diameter

area_inner = np.pi*(d_inner/2)**2
area_outer = np.pi*(d_outer/2)**2 - np.pi*(d_dielectric/2)**2
# Represent the annular dielectric as an equivalent rectangular slab (see below)
t_dielectric = (d_dielectric - d_inner)/2
w_dielectric = t_dielectric * 2*np.pi/np.log(d_dielectric/d_inner)

result = solve_coupled_conductor_heat_flow(
    k_inner=k_ss,                 # inner conductor material
    k_outer=k_ss,                 # outer conductor material
    k_dielectric=k_ptfe,          # dielectric material
    area_inner=area_inner,        # m^2
    area_outer=area_outer,        # m^2
    length=0.1,                   # m
    dielectric_thickness=t_dielectric,  # m
    dielectric_width=w_dielectric,      # m
    inner_anchors=[(0, 300), (0.1, 4)],            # (position, temperature)
    outer_anchors=[(0, 300), (0.04, 40), (0.1, 4)], # outer conductor heat-sunk at 3 stages
)

for T, P in sorted(result['stage_powers'].items()):
    print(f'{T:6.1f} K stage: {P*1e3:+8.4f} mW deposited')

giving us

   4.0 K stage:  +2.5506 mW deposited
  40.0 K stage: +119.1332 mW deposited
 300.0 K stage: -121.6839 mW deposited

The result dictionary also contains the grid x, the temperature profiles T_inner and T_outer, and per-anchor powers in anchor_powers. Positive power means heat flows from the cable into that stage; the stage powers always sum to ~0 (energy conservation, reported as power_balance_error).

Model: PDE, coupling, and boundary conditions

The inner and outer conductors are each treated as a 1D rod with temperature-dependent conductivity, coupled to each other through the dielectric at every point along the length. In steady state:

d/dx [ k_i(T_i) A_i dT_i/dx ] + g(T_i, T_o) (T_o - T_i) = 0     (inner conductor)
d/dx [ k_o(T_o) A_o dT_o/dx ] + g(T_i, T_o) (T_i - T_o) = 0     (outer conductor)

where T_i(x) and T_o(x) are the conductor temperatures, A_i and A_o are their cross-sectional areas, and g is the dielectric coupling conductance per unit length (W/m/K):

g(T_i, T_o) = (w/t) * k_d_eff,   with   k_d_eff = (1/(T_o - T_i)) ∫ k_d(T) dT  from T_i to T_o

The dielectric's own temperature profile is not solved for — it is treated as a pure thermal link between the two conductors (its axial conduction is negligible compared to the metal conductors, and its radial heat capacity doesn't matter in steady state).

Boundary conditions: Every anchor is a Dirichlet (fixed-temperature) condition, T(x_a) = T_a, applied at the grid node nearest the requested position. Any conductor end that has no anchor is naturally adiabatic (zero heat flux). The power reported for each anchor is the net conducted heat arriving at that node from its neighbors and from the opposite conductor through the dielectric.

Rectangular-slab dielectric approximation: The solver models the dielectric as a rectangular slab of width w and thickness t (heat-flow path length between the conductors). To represent a true annular dielectric with inner diameter a and outer diameter b, use the exact per-unit-length equivalence of an annulus, w/t = 2π/ln(b/a) — e.g. t = (b-a)/2 and w = t*2π/ln(b/a) as in the example above. Any (w, t) pair with the correct ratio gives identical results.

Numerics: The equations are discretized with finite volumes on a fixed uniform grid (npoints nodes, no adaptive meshing). The conductance of each grid link is computed from the thermal conductivity integral between the two node temperatures, so the strong temperature dependence of k(T) is captured exactly along each link rather than sampled at one point (the scheme is 2nd-order accurate: doubling npoints reduces profile errors ~4x). The nonlinear system is solved by Picard iteration: conductances are frozen, the resulting linear system is solved directly (sparse LU), conductances are re-evaluated at the new temperatures, and the loop repeats until the largest temperature change is below tol. Because every linear system is a pure resistor network with positive conductances and no sources, the discrete maximum principle guarantees each iterate lies between the coldest and hottest anchor temperatures — temperatures can never go negative or overshoot the boundary conditions, and k(T) is never evaluated outside the anchored temperature range. If convergence stalls for a very sharply peaked k(T), set relaxation < 1 (e.g. 0.7).

Verifying the solver

verify_coax.py in the repository root checks the solver against simpler known solutions, which are also good sanity checks for your own configurations:

  1. Constant k, zero coupling (k_dielectric = 0): each conductor reduces to an independent rod — profiles must be exactly linear and the power must equal k*A*(T1-T2)/L.
  2. Nonlinear k, zero coupling: the end powers must match the thermal conductivity integral result from calculate_thermal_transfer().
  3. Fin equation: with constant conductivities and the outer conductor clamped isothermal at T0 (anchored at every node), the inner conductor obeys the classic fin equation d²T/dx² = m²(T - T0) with m² = (k_d w/t)/(k_i A_i), which has an analytic sinh/cosh solution for both the profile and the end heat flows.
  4. Conservation and convergence: for any configuration, the stage powers must sum to zero, and the answers must stop changing as npoints is increased (refine until the quantity you care about changes by less than your tolerance).

Plotting Thermal Conductivity Curves

image

Plotting code here: https://github.com/amccaugh/cryoheatflow/blob/main/plot_thermal_conductivities.py

Available Materials

Thermal Conductivity Functions

The package provides thermal conductivity functions for various materials. Most are sourced from the NIST cryogenic thermal conductivity reference; tabulated data materials are noted separately.

  • k_ss - Stainless steel (316/314/304L), valid 1–300 K
  • k_cuni - 70-30 CuNi cupronickel, valid 1–300 K
  • k_al6061 - Aluminum 6061-T6, valid 1–300 K
  • k_al6063 - Aluminum 6063-T5, valid 4–300 K
  • k_al1100 - Aluminum 1100, valid 4–300 K
  • k_brass - Brass (UNS C26000), valid 5–110 K
  • k_becu - Beryllium copper, valid 4–120 K
  • k_cu_rrr50 - Copper (RRR=50, typically ETP or OFHC), valid 4–300 K
  • k_cu_rrr100 - Copper (RRR=100), valid 4–300 K
  • k_g10 - Fiberglass-epoxy (G-10), valid 4–300 K
  • k_nylon - Nylon (polyamide), valid 4–300 K
  • k_ptfe - PTFE (Teflon), valid 4–300 K
  • k_phosphor_bronze - Phosphor bronze (94.8% Cu, 5% Sn, 0.2% P) — tabulated data from Lake Shore Cryotronics, valid 1–300 K
  • k_nichrome - Nichrome (80% Ni, 20% Cr) — tabulated data from Lake Shore Cryotronics, valid 4–300 K
  • k_manganin - Manganin (83% Cu, 13% Mn, 4% Ni) — tabulated data from Lake Shore Cryotronics, valid 0.1–300 K

Each function's valid range is also available programmatically as k_fun.T_min / k_fun.T_max.

Extrapolating Beyond the Data Range

Each conductivity function returns NaN outside its measured data range (listed above). If you need values slightly beyond that range — say a BeCu part of a cable that runs up to 300K when the BeCu data stops at 120K — you can opt in to an artificial power-law extension with extrapolate():

from cryoheatflow import k_becu, extrapolate

# BeCu data spans 4-120 K; stretch it to cover 1-300 K
k = extrapolate(k_becu, T_min=1, T_max=300)

k(200)   # power-law continuation of the curve
# ExtrapolationWarning: k_becu evaluated outside its 4-120 K data range; ...

k(0.5)   # still NaN -- beyond the limit you stated

The extension continues the curve's log-log slope at the data boundary as a power law, which is smooth and physically motivated (many materials follow approximate power laws at cryogenic temperatures). An ExtrapolationWarning is issued whenever extrapolated values are actually used, and temperatures beyond your stated T_min/T_max still return NaN. The wrapped function works anywhere the original does, including solve_coupled_conductor_heat_flow() and calculate_thermal_transfer().

Caution: the extension is a guess, not data — especially below a few kelvin, where real materials can deviate from a single power law by large factors.

Thermal Boundary Conductance Functions

  • h_grease - Thermal conductance of grease for given contact area
  • h_solder_pb_sn - Thermal conductance of standard lead-tin (PbSn) solder for given contact area

Emissivity Values

  • Al_polished - Polished aluminum (ε = 0.03)
  • Al_oxidized - Oxidized aluminum (ε = 0.3)
  • Cu_polished - Polished copper (ε = 0.02)
  • Cu_oxidized - Oxidized copper (ε = 0.6)
  • brass_polished - Polished brass (ε = 0.03)
  • brass_oxidized - Oxidized brass (ε = 0.6)
  • stainless - Stainless steel (ε = 0.07)
  • mylar - Mylar (ε = 0.05)

Area Calculations

The package includes functions for calculating cross-sectional areas:

  • tube_area(diameter, wall_thickness) - Annular cross-section area
  • cylinder_area(diameter) - Circular cross-section area
  • wire_gauge_area(awg) - Wire cross-section area based on AWG (American Wire Gauge)
  • wire_swg_area(swg) - Wire cross-section area based on SWG (Standard/Imperial Wire Gauge, BS 3737), valid for gauges 1–50
  • coax_141, coax_085, coax_047, coax_034 - Predefined coaxial cable areas

Data Sources

Thermal conductivity data is sourced from:

Emissivity and thermal boundary conductance values are from Ekin, J. (2006), Experimental Techniques for Low-Temperature Measurements, Oxford University Press, Oxford, UK.

Requirements

  • Python >= 3
  • NumPy
  • SciPy
  • Matplotlib (for plotting examples)

Acknowledgements

This package was developed by Adam McCaughan. Special thanks to the wider cryogenic-science community for the invaluable data used in this package. If you use this package in your work, please consider citing the relevant sources and acknowledging the authors.

About

Cryogenic thermal calculations for cryostat design

Resources

Stars

8 stars

Watchers

1 watching

Forks

Releases

Packages

Contributors

Languages