The specific heat of a gas mixture¶

The specific heat of a gas mixture, the mass-weighted mean of its gases', c_p = sum w_i c_p,i with w_i = x_i M_i / M the mass fractions, x_i = n_i / sum n the amount fractions and M the mixture's molar mass (lusine.properties.molar_mass's, passed in): an ideal mixture, each gas holding its heat as it would alone.

Over a game's whole table of parts, as molar_mass is: each part's mixture, its amount (any amounts: only their proportions count) and its gas, a relation to whatever table the game's column relates to, whose molar_mass and specific_heat it reads. It folds the parts of the item's own mixture; each sum is added in the fixed tree. Source: Cengel & Boles, Thermodynamics: An Engineering Approach, ch. 13 (gas mixtures). Checked against a reference in numpy, which has its own tests.

Source

Cengel & Boles, Thermodynamics, ch. 13 (gas mixtures).

Proof

reference, limit, invariant.

mixture.trait

namespace lusine.properties;

float<J/(kg K)> mixture(list<struct { relation mixture; float fraction; relation gas; }> parts,
                        float<kg/mol> molar_mass)
{
  let total = sum(fraction) where mixture == it;
  let heat = sum(fraction * gas.molar_mass * gas.specific_heat) where mixture == it;
  return heat / (total * molar_mass);
}

A gas mixture's c_p, the mass-weighted mean of its gases', c_p = Σ w_i c_p,i with w_i = x_i M_i / M. The kernel takes a game's whole table of parts, list<struct { relation mixture; float fraction; relation gas; }>, as lusine.properties.molar_mass does, and the mixture's molar mass, which that kernel gives.

In [1]:

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'))
sys.path.insert(0, os.path.join(root, 'kernels'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
from lusine import rows
from lusine.kernels.properties._gases import DRY_AIR, GASES
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24

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')) sys.path.insert(0, os.path.join(root, 'kernels')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run from lusine import rows from lusine.kernels.properties._gases import DRY_AIR, GASES F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24

The kernel¶

The law as the library states it: the fields it needs of a part, and of a part's gas its molar mass and specific heat.

In [2]:

print(trait_run.source('properties/mixture.trait'))

print(trait_run.source('properties/mixture.trait'))

// lusine.properties.mixture: the specific heat of a gas mixture, the mass-weighted mean of its gases', c_p = sum w_i
// c_p,i with w_i = x_i M_i / M the mass fractions, x_i = n_i / sum n the amount fractions and M the mixture's molar mass
// (lusine.properties.molar_mass's, passed in): an ideal mixture, each gas holding its heat as it would alone. Over a
// game's whole table of parts, as molar_mass is: each part's mixture, its amount (any amounts: only
// their proportions count) and its gas, a relation to whatever table the game's column relates to, whose molar_mass
// and specific_heat it reads. It folds the parts of the item's own mixture; each sum is added in the fixed tree.
// Source: Cengel & Boles, Thermodynamics: An Engineering Approach, ch. 13 (gas mixtures).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.properties;

float<J/(kg K)> mixture(list<struct { relation mixture; float fraction; relation gas; }> parts,
                        float<kg/mol> molar_mass)
{
  let total = sum(fraction) where mixture == it;
  let heat = sum(fraction * gas.molar_mass * gas.specific_heat) where mixture == it;
  return heat / (total * molar_mass);
}

The program¶

A game's tables, written out: its gases (each a molar mass and a c_p near 300 K, the CRC Handbook's and Engineering ToolBox's, as the reference's), its mixtures and its parts, each part a mixture's amount of a gas: dry air, the same air in thousands of moles, nitrogen alone, equal amounts of nitrogen and steam, and air with its water rising from none to 5% of the moles. main() binds molar_mass over the mixtures, then this kernel with it; the planted fault below binds the page's own wrong kernel in its place.

In [3]:

WATER = [round(0.005 * i, 3) for i in range(11)]
MIXTURES = {'air': DRY_AIR, 'kiloair': {g: 1000 * x for g, x in DRY_AIR.items()}, 'nitrogen': {'N2': 1.0},
            'even': {'N2': 1.0, 'H2O': 1.0}}
MIXTURES.update({f'humid{i}': dict(DRY_AIR, H2O=w / (1 - w)) for i, w in enumerate(WATER)})


def program(kernel='lusine.properties.mixture', extra=''):
    gases = ',\n'.join(f'  {{ name = "{g}", molar_mass = {r["molar_mass"]!r} kg/mol, '
                       f'specific_heat = {r["c_p"]!r} J/(kg K) }}' for g, r in GASES.items())
    mixtures = ', '.join(f'{{ name = "{m}" }}' for m in MIXTURES)
    parts = ',\n'.join(f'  {{ mixture = mixtures["{m}"], fraction = {float(x)!r}, gas = gases["{g}"] }}'
                       for m, c in MIXTURES.items() for g, x in c.items())
    return f'''namespace mix;

lusine.version version = {{ {{ number = 1 }} }};

struct gas {{ text name; float<kg/mol> molar_mass; float<J/(kg K)> specific_heat; }};
table<gas> gases(capacity = 8, key = name) = {{
{gases},
}};
struct mixture {{ text name; float<kg/mol> molar_mass; float<J/(kg K)> specific_heat; }};
table<mixture> mixtures(capacity = 16, key = name) = {{ {mixtures} }};
struct part {{ relation<mixtures> mixture; float fraction; relation<gases> gas; }};
table<part> parts(capacity = 96) = {{
{parts},
}};
{extra}
void main()
{{
  mixtures.molar_mass = lusine.properties.molar_mass(parts);
  mixtures.specific_heat = {kernel}(parts, mixtures.molar_mass);
}}
'''

LAW = program()
print(LAW)


def heats(text=LAW):
    """{mixture: (its molar mass, its c_p)}, and the run."""
    run = trait_run.Program(['properties/mixture.trait', 'properties/molar_mass.trait'], text, {})
    m, c = run.steps('mix.main', 1, [('mix.mixtures', 'molar_mass'), ('mix.mixtures', 'specific_heat')])
    return dict(zip(MIXTURES, zip(m[0], c[0]))), run

C, run = heats()

WATER = [round(0.005 * i, 3) for i in range(11)] MIXTURES = {'air': DRY_AIR, 'kiloair': {g: 1000 * x for g, x in DRY_AIR.items()}, 'nitrogen': {'N2': 1.0}, 'even': {'N2': 1.0, 'H2O': 1.0}} MIXTURES.update({f'humid{i}': dict(DRY_AIR, H2O=w / (1 - w)) for i, w in enumerate(WATER)}) def program(kernel='lusine.properties.mixture', extra=''): gases = ',\n'.join(f' {{ name = "{g}", molar_mass = {r["molar_mass"]!r} kg/mol, ' f'specific_heat = {r["c_p"]!r} J/(kg K) }}' for g, r in GASES.items()) mixtures = ', '.join(f'{{ name = "{m}" }}' for m in MIXTURES) parts = ',\n'.join(f' {{ mixture = mixtures["{m}"], fraction = {float(x)!r}, gas = gases["{g}"] }}' for m, c in MIXTURES.items() for g, x in c.items()) return f'''namespace mix; lusine.version version = {{ {{ number = 1 }} }}; struct gas {{ text name; float<kg/mol> molar_mass; float<J/(kg K)> specific_heat; }}; table<gas> gases(capacity = 8, key = name) = {{ {gases}, }}; struct mixture {{ text name; float<kg/mol> molar_mass; float<J/(kg K)> specific_heat; }}; table<mixture> mixtures(capacity = 16, key = name) = {{ {mixtures} }}; struct part {{ relation<mixtures> mixture; float fraction; relation<gases> gas; }}; table<part> parts(capacity = 96) = {{ {parts}, }}; {extra} void main() {{ mixtures.molar_mass = lusine.properties.molar_mass(parts); mixtures.specific_heat = {kernel}(parts, mixtures.molar_mass); }} ''' LAW = program() print(LAW) def heats(text=LAW): """{mixture: (its molar mass, its c_p)}, and the run.""" run = trait_run.Program(['properties/mixture.trait', 'properties/molar_mass.trait'], text, {}) m, c = run.steps('mix.main', 1, [('mix.mixtures', 'molar_mass'), ('mix.mixtures', 'specific_heat')]) return dict(zip(MIXTURES, zip(m[0], c[0]))), run C, run = heats()

namespace mix;

lusine.version version = { { number = 1 } };

struct gas { text name; float<kg/mol> molar_mass; float<J/(kg K)> specific_heat; };
table<gas> gases(capacity = 8, key = name) = {
  { name = "N2", molar_mass = 0.0280134 kg/mol, specific_heat = 1040.0 J/(kg K) },
  { name = "O2", molar_mass = 0.0319988 kg/mol, specific_heat = 918.0 J/(kg K) },
  { name = "Ar", molar_mass = 0.039948 kg/mol, specific_heat = 520.3 J/(kg K) },
  { name = "CO2", molar_mass = 0.0440095 kg/mol, specific_heat = 844.0 J/(kg K) },
  { name = "H2O", molar_mass = 0.01801528 kg/mol, specific_heat = 1864.0 J/(kg K) },
};
struct mixture { text name; float<kg/mol> molar_mass; float<J/(kg K)> specific_heat; };
table<mixture> mixtures(capacity = 16, key = name) = { { name = "air" }, { name = "kiloair" }, { name = "nitrogen" }, { name = "even" }, { name = "humid0" }, { name = "humid1" }, { name = "humid2" }, { name = "humid3" }, { name = "humid4" }, { name = "humid5" }, { name = "humid6" }, { name = "humid7" }, { name = "humid8" }, { name = "humid9" }, { name = "humid10" } };
struct part { relation<mixtures> mixture; float fraction; relation<gases> gas; };
table<part> parts(capacity = 96) = {
  { mixture = mixtures["air"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["air"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["air"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["air"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["kiloair"], fraction = 780.84, gas = gases["N2"] },
  { mixture = mixtures["kiloair"], fraction = 209.46, gas = gases["O2"] },
  { mixture = mixtures["kiloair"], fraction = 9.34, gas = gases["Ar"] },
  { mixture = mixtures["kiloair"], fraction = 0.412, gas = gases["CO2"] },
  { mixture = mixtures["nitrogen"], fraction = 1.0, gas = gases["N2"] },
  { mixture = mixtures["even"], fraction = 1.0, gas = gases["N2"] },
  { mixture = mixtures["even"], fraction = 1.0, gas = gases["H2O"] },
  { mixture = mixtures["humid0"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid0"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid0"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid0"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid0"], fraction = 0.0, gas = gases["H2O"] },
  { mixture = mixtures["humid1"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid1"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid1"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid1"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid1"], fraction = 0.005025125628140704, gas = gases["H2O"] },
  { mixture = mixtures["humid2"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid2"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid2"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid2"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid2"], fraction = 0.010101010101010102, gas = gases["H2O"] },
  { mixture = mixtures["humid3"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid3"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid3"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid3"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid3"], fraction = 0.015228426395939085, gas = gases["H2O"] },
  { mixture = mixtures["humid4"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid4"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid4"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid4"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid4"], fraction = 0.020408163265306124, gas = gases["H2O"] },
  { mixture = mixtures["humid5"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid5"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid5"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid5"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid5"], fraction = 0.025641025641025644, gas = gases["H2O"] },
  { mixture = mixtures["humid6"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid6"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid6"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid6"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid6"], fraction = 0.030927835051546393, gas = gases["H2O"] },
  { mixture = mixtures["humid7"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid7"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid7"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid7"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid7"], fraction = 0.036269430051813475, gas = gases["H2O"] },
  { mixture = mixtures["humid8"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid8"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid8"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid8"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid8"], fraction = 0.04166666666666667, gas = gases["H2O"] },
  { mixture = mixtures["humid9"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid9"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid9"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid9"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid9"], fraction = 0.04712041884816754, gas = gases["H2O"] },
  { mixture = mixtures["humid10"], fraction = 0.78084, gas = gases["N2"] },
  { mixture = mixtures["humid10"], fraction = 0.20946, gas = gases["O2"] },
  { mixture = mixtures["humid10"], fraction = 0.00934, gas = gases["Ar"] },
  { mixture = mixtures["humid10"], fraction = 0.000412, gas = gases["CO2"] },
  { mixture = mixtures["humid10"], fraction = 0.052631578947368425, gas = gases["H2O"] },
};

void main()
{
  mixtures.molar_mass = lusine.properties.molar_mass(parts);
  mixtures.specific_heat = lusine.properties.mixture(parts, mixtures.molar_mass);
}

Reference: dry air¶

The International Standard Atmosphere's air has c_p near 1005 J/(kg K); from these gases' values the float64 reference gives 1004.96, and the kernel the same within binary32's rounding.

In [4]:

from lusine.kernels.properties.molar_mass import molar_mass as reference_m
from lusine.kernels.properties.mixture import mixture as reference
cp = C['air'][1]
print(f'{cp:.2f} J/(kg K); the float64 reference {reference(DRY_AIR, reference_m(DRY_AIR)):.2f}')
assert abs(cp - 1005.0) < 10.0

from lusine.kernels.properties.molar_mass import molar_mass as reference_m from lusine.kernels.properties.mixture import mixture as reference cp = C['air'][1] print(f'{cp:.2f} J/(kg K); the float64 reference {reference(DRY_AIR, reference_m(DRY_AIR)):.2f}') assert abs(cp - 1005.0) < 10.0

1004.96 J/(kg K); the float64 reference 1004.96

Reference: bit for bit, the folds in binary32¶

Every mixture's c_p is a binary32 model's bit for bit: the mixture's parts in slot order, their amounts and their amount × M_i × c_p,i each added in Trait's fixed tree, and the second divided by the first times the mixture's molar mass (the program's, held bit for bit on the molar mass's own page). The float64 reference agrees within 12 × 2⁻²⁴.

In [5]:

def tree(values):
    if not values:
        return F32(0.0)
    if len(values) == 1:
        return F32(values[0])
    mid = (len(values) + 1) // 2
    return F32(tree(values[:mid]) + tree(values[mid:]))

cells = {r['key']: r for r in rows.rows(run.unit, 'mix.gases')}


def model(composition, molar_mass):
    total = tree([F32(x) for x in composition.values()])
    heat = tree([F32(F32(F32(x) * F32(cells[g]['molar_mass'])) * F32(cells[g]['specific_heat']))
                 for g, x in composition.items()])
    return F32(heat / F32(total * F32(molar_mass)))


def model_holds(C):
    return all(F32(C[m][1]) == model(c, C[m][0]) for m, c in MIXTURES.items())

assert model_holds(C)
worst = max(abs(C[m][1] - reference(c, reference_m(c))) / C[m][1] for m, c in MIXTURES.items())
assert worst <= 12 * U, worst / U
print(f'{len(MIXTURES)} mixtures bit for bit with the binary32 model; the float64 reference within {worst / U:.1f} x 2^-24')

def tree(values): if not values: return F32(0.0) if len(values) == 1: return F32(values[0]) mid = (len(values) + 1) // 2 return F32(tree(values[:mid]) + tree(values[mid:])) cells = {r['key']: r for r in rows.rows(run.unit, 'mix.gases')} def model(composition, molar_mass): total = tree([F32(x) for x in composition.values()]) heat = tree([F32(F32(F32(x) * F32(cells[g]['molar_mass'])) * F32(cells[g]['specific_heat'])) for g, x in composition.items()]) return F32(heat / F32(total * F32(molar_mass))) def model_holds(C): return all(F32(C[m][1]) == model(c, C[m][0]) for m, c in MIXTURES.items()) assert model_holds(C) worst = max(abs(C[m][1] - reference(c, reference_m(c))) / C[m][1] for m, c in MIXTURES.items()) assert worst <= 12 * U, worst / U print(f'{len(MIXTURES)} mixtures bit for bit with the binary32 model; the float64 reference within {worst / U:.1f} x 2^-24')

15 mixtures bit for bit with the binary32 model; the float64 reference within 1.4 x 2^-24

Limit: one gas alone¶

A mixture of one gas has that gas's c_p, within an ulp (fl(M c_p) / M).

In [6]:

n2 = F32(GASES['N2']['c_p'])
assert abs(C['nitrogen'][1] - n2) <= 2 * U * n2
print('nitrogen alone:', C['nitrogen'][1])

n2 = F32(GASES['N2']['c_p']) assert abs(C['nitrogen'][1] - n2) <= 2 * U * n2 print('nitrogen alone:', C['nitrogen'][1])

nitrogen alone: 1040.0

Invariant: by mass, not by amount¶

Equal amounts of nitrogen and steam weigh each c_p by its gas's mass, not by its amount; the same air in thousands of moles holds the same heat; and humid air holds more heat per kilogram as its water rises.

In [7]:

N2, H2O = GASES['N2'], GASES['H2O']
by_mass = (N2['molar_mass'] * N2['c_p'] + H2O['molar_mass'] * H2O['c_p']) / (N2['molar_mass'] + H2O['molar_mass'])


def mass_holds(C):
    return abs(C['even'][1] - by_mass) <= 8 * U * by_mass and abs(C['kiloair'][1] - C['air'][1]) <= 8 * U * C['air'][1]

humid = np.array([C[f'humid{i}'][1] for i in range(len(WATER))])
assert mass_holds(C) and np.all(np.diff(humid) > 0)
print(f'even nitrogen and steam: {C["even"][1]:.2f} J/(kg K), by mass {by_mass:.2f}; by amount it would be '
      f'{(N2["c_p"] + H2O["c_p"]) / 2:.0f}')
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(WATER, humid, 'o-')
ax.set_xlabel('water, amount fraction')
ax.set_ylabel('c_p (J/(kg K))')
plt.show()

N2, H2O = GASES['N2'], GASES['H2O'] by_mass = (N2['molar_mass'] * N2['c_p'] + H2O['molar_mass'] * H2O['c_p']) / (N2['molar_mass'] + H2O['molar_mass']) def mass_holds(C): return abs(C['even'][1] - by_mass) <= 8 * U * by_mass and abs(C['kiloair'][1] - C['air'][1]) <= 8 * U * C['air'][1] humid = np.array([C[f'humid{i}'][1] for i in range(len(WATER))]) assert mass_holds(C) and np.all(np.diff(humid) > 0) print(f'even nitrogen and steam: {C["even"][1]:.2f} J/(kg K), by mass {by_mass:.2f}; by amount it would be ' f'{(N2["c_p"] + H2O["c_p"]) / 2:.0f}') fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(WATER, humid, 'o-') ax.set_xlabel('water, amount fraction') ax.set_ylabel('c_p (J/(kg K))') plt.show()

even nitrogen and steam: 1362.51 J/(kg K), by mass 1362.51; by amount it would be 1452

Planted fault¶

A wrong kernel over the same list, written here in Trait and bound in the library's place: the amount-weighted mean of the gases' c_p, the mass fractions forgotten. The mass check and the binary32 model each catch it.

In [8]:

BY_AMOUNT = program('by_amount', '''
float<J/(kg K)> by_amount(list<struct { relation mixture; float fraction; relation gas; }> parts, float<kg/mol> molar_mass)
{
  let total = sum(fraction) where mixture == it;
  return (sum(fraction * gas.specific_heat) where mixture == it) / total;
}
''')
W, _ = heats(BY_AMOUNT)
assert not mass_holds(W) and not model_holds(W)
print('PLANTED: the amount-weighted mean is caught')

BY_AMOUNT = program('by_amount', ''' float<J/(kg K)> by_amount(list<struct { relation mixture; float fraction; relation gas; }> parts, float<kg/mol> molar_mass) { let total = sum(fraction) where mixture == it; return (sum(fraction * gas.specific_heat) where mixture == it) / total; } ''') W, _ = heats(BY_AMOUNT) assert not mass_holds(W) and not model_holds(W) print('PLANTED: the amount-weighted mean is caught')

PLANTED: the amount-weighted mean is caught