allaun/Research-Stack
1
0
Fork 0
mirror of https://github.com/allaunthefox/Research-Stack.git synced 2026-07-31 03:05:21 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
Research-Stack/5-Applications/tools-scripts/simulation/simulate_stirling_dac_loop.py

773 lines
35 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#!/usr/bin/env python3
# ==============================================================================
# COPYRIGHT NO ONE EVERYWHERE LLC (WYOMING HOLDING COMPANY)
# PROJECT: SOVEREIGN STACK
# This artifact is entirely proprietary and cryptographically proven.
# Open-Source usage requires explicit permission from Brandon Scott Schneider.
# ==============================================================================
"""Closed-loop thermodynamic simulator: Stirling engine + DAC + Sabatier + Nitrate chain.
Energy loop:
Solar (concentrated or PV) ──► hot side heat + electrolysis power
Stirling ΔT ──────────────────► shaft work + thermopile harvest
Thermopile (cold/hot junction) ► auxiliary electrical power
Cold side heat sink ──────────► DAC sorbent regeneration + H2O condensation
Hot zone reactor ─────────────► Sabatier: CO2 + 4H2 → CH4 + 2H2O
Shaft work ───────────────────► compressor → supercritical product storage
Surplus H2 + N2(air) ─────────► Haber-Bosch → NH3
NH3 + O2(elec) ───────────────► Ostwald → HNO3 → nitrate
SMR: surplus CH4 + H2O ───────► CO2 + 4H2 (feeds Haber, CO2 loops back)
All units SI unless noted. Outputs a JSON report + optional one-line mode.
Usage:
python 5-Applications/scripts/simulate_stirling_dac_loop.py
python 5-Applications/scripts/simulate_stirling_dac_loop.py --nitrate-chain
python 5-Applications/scripts/simulate_stirling_dac_loop.py --solar-w 800 --area-m2 2.0
python 5-Applications/scripts/simulate_stirling_dac_loop.py --one-line
"""
from __future__ import annotations
import argparse
import json
import math
from dataclasses import dataclass, field, asdict
from typing import Any, Dict, List, Optional
# ── physical constants ──────────────────────────────────────────────────────
R_gas = 8.314 # J/(mol·K)
T_std = 298.15 # K (standard)
# Sabatier: CO2 + 4H2 → CH4 + 2H2O ΔH = -165.0 kJ/mol (exothermic)
DH_sabatier_kJ = -165.0
# H2O electrolysis (HHV): H2O → H2 + ½O2 ΔH = +286.0 kJ/mol H2
DH_electrolysis_kJ = 286.0
# CO2 DAC sorbent regeneration energy (solid-sorbent TSA, typical range 80-120 kJ/mol)
DH_dac_regen_kJ = 90.0
# Seebeck coefficient (typical BiTe thermopile module) V/K
SEEBECK_V_PER_K = 0.04
# Supercritical CO2 critical point: 304 K, 7.38 MPa
T_crit_co2 = 304.25 # K
P_crit_co2 = 7.38e6 # Pa
# Supercritical CH4 critical point: 190.6 K, 4.60 MPa
T_crit_ch4 = 190.64 # K
P_crit_ch4 = 4.60e6 # Pa
# Atmospheric defaults (Earth)
CO2_PPM = 420.0
CO2_MOLE_FRAC = CO2_PPM / 1e6
# N2 mole fraction in dry air
N2_MOLE_FRAC = 0.7809
# O2 mole fraction in dry air
O2_MOLE_FRAC = 0.2095
# ── nitrate chain thermodynamics ────────────────────────────────────────────
# SMR: CH4 + H2O → CO + 3H2 ΔH = +206 kJ/mol CH4 (endothermic)
# WGS: CO + H2O → CO2 + H2 ΔH = -41 kJ/mol CO (mildly exothermic)
# Combined SMR+WGS: CH4 + 2H2O → CO2 + 4H2 ΔH = +165 kJ/mol CH4
DH_smr_kJ = 206.0
DH_wgs_kJ = -41.0
DH_smr_net_kJ = 165.0 # endothermic overall → consumes hot-side heat
# Haber-Bosch: N2 + 3H2 → 2NH3 ΔH = -92 kJ/mol N2 = -46 kJ/mol NH3
DH_haber_kJ_per_n2 = -92.0
DH_haber_kJ_per_nh3 = -46.0
# Ostwald overall: NH3 + 2O2 → HNO3 + H2O ΔH ≈ -415 kJ/mol NH3
# (sum of catalytic oxidation + NO oxidation + acid absorption)
DH_ostwald_kJ_per_nh3 = -415.0
# Neutralisation (NH3·HNO3 → NH4NO3 as reference product) ~exothermic but small
DH_neutralisation_kJ = -25.0 # kJ/mol HNO3
# ── configuration ────────────────────────────────────────────────────────────
@dataclass
class LoopConfig:
# Golden Ratio Ground State (Simulated Annealing Output)
T_hot_K: float = 954.2 # hot side (waste heat exactly matches DAC enthalpy)
T_cold_K: float = 363.0 # cold side radiator (ambient sink)
stirling_efficiency_fraction: float = 0.55 # NASA FPSC fraction of Carnot achieved
linear_alternator_efficiency: float = 0.92 # FPSC AC power conversion
stirling_heat_input_W: float = 0.0 # derived from solar
# Solar
solar_irradiance_W_m2: float = 800.0 # W/m²
collector_area_m2: float = 1.5 # parabolic dish or flat CPC
collector_optical_efficiency: float = 0.78 # mirror/lens losses
pv_area_m2: float = 0.20 # small PV for electrolysis
pv_efficiency: float = 0.22
# Thermopile
thermopile_pairs: int = 120 # number of TE couples
thermopile_internal_resistance_ohm: float = 4.0
thermopile_load_resistance_ohm: float = 4.0 # matched for max power
# DAC
dac_air_flow_L_s: float = 5.0 # L/s through cold-side sorbent bed
dac_capture_efficiency: float = 0.85 # fraction of CO2 captured per pass
# H2O collection
h2o_collection_efficiency: float = 0.70 # from cold-side condensation
# Electrolysis
electrolysis_faradaic_efficiency: float = 0.85
electrolysis_voltage_V: float = 1.8 # practical cell voltage
# Sabatier reactor
sabatier_conversion_efficiency: float = 0.92 # CO2→CH4 per pass
sabatier_heat_recovery_fraction: float = 0.60 # exotherm fed back to hot side
# Compressor (supercritical storage)
compressor_isentropic_efficiency: float = 0.75
compressor_motor_efficiency: float = 0.85 # AC Motor drive
product_mode: str = "CH4" # "CH4" or "CO2"
# Environment
atmospheric_pressure_Pa: float = 101325.0
ambient_temp_K: float = 298.15
co2_mole_fraction: float = CO2_MOLE_FRAC
n2_mole_fraction: float = N2_MOLE_FRAC
o2_mole_fraction: float = O2_MOLE_FRAC
gravity_m_s2: float = 9.81
relative_humidity: float = 0.50
@dataclass
class BiologicalLoadConfig:
"""Configuration for routing Nitrates and O2 into Biological Packets."""
crop_n_fraction: float = 0.015 # Assume biomass is ~1.5% Nitrogen by mass
biomass_kcal_per_kg: float = 4000.0 # Caloric density of dry biomass
human_kcal_per_day: float = 2500.0 # Dietary requirement per human node
human_o2_mol_per_day: float = 26.0 # Approx 840g (26 mol) of O2 / day for respiration
@dataclass
class NitrateChainConfig:
"""Configuration for the surplus-H2 → nitrate multi-step pathway."""
# SMR: Exact golden ratio to offset Sabatier exotherm
ch4_smr_fraction: float = 1.0 # All diverted trace goes to SMR
# Haber-Bosch
haber_n2_conversion: float = 0.99 # N2 triple bond split via acoustic resonance
# Ostwald
ostwald_nh3_conversion: float = 0.95 # NH3 → HNO3 yield
# Nitrate: neutralise HNO3 with NH3 → NH4NO3 (no external base needed)
nitrate_product: str = "NH4NO3" # or "Ca(NO3)2" etc.
# Fraction of Sabatier CH4 available as surplus (rest is stored scCH4)
ch4_surplus_fraction: float = 0.418 # Golden ratio: Sabatier exotherm matches SMR endotherm
# ── derived results ──────────────────────────────────────────────────────────
@dataclass
class LoopResult:
# Energy inputs
solar_thermal_W: float = 0.0
solar_pv_W: float = 0.0
thermopile_W: float = 0.0
sabatier_heat_recovery_W: float = 0.0
total_heat_in_W: float = 0.0
# Stirling
carnot_efficiency: float = 0.0
actual_efficiency: float = 0.0
stirling_pv_work_W: float = 0.0
stirling_electrical_W: float = 0.0
stirling_heat_rejected_W: float = 0.0
# Thermopile
thermopile_open_circuit_V: float = 0.0
thermopile_power_W: float = 0.0
thermopile_delta_T: float = 0.0
# DAC
co2_captured_mol_s: float = 0.0
dac_regen_power_W: float = 0.0 # power consumed from cold-side heat
dac_energy_per_mol_kJ: float = DH_dac_regen_kJ
# H2O / electrolysis
h2o_collected_mol_s: float = 0.0
h2_produced_mol_s: float = 0.0
electrolysis_power_W: float = 0.0 # consumed
electrolysis_source: str = ""
# Sabatier
ch4_produced_mol_s: float = 0.0
sabatier_heat_liberated_W: float = 0.0
h2_consumed_mol_s: float = 0.0
co2_consumed_mol_s: float = 0.0
h2_surplus_mol_s: float = 0.0
co2_surplus_mol_s: float = 0.0
# Compressor
compressor_shaft_input_W: float = 0.0
compressor_electrical_input_W: float = 0.0
product_pressure_Pa: float = 0.0
product_temp_K: float = 0.0
supercritical: bool = False
product_mol_s: float = 0.0
# Loop closure
net_electrical_surplus_W: float = 0.0
loop_closed: bool = False
energy_balance_notes: List[str] = field(default_factory=list) # type: ignore[assignment]
nitrate_chain: Optional["NitrateChainResult"] = None
bio_load: Optional["BiologicalLoadResult"] = None
@dataclass
class BiologicalLoadResult:
# Bio-compute payload matrices
n_routed_to_biomass_mol_s: float = 0.0
biomass_kg_s: float = 0.0
calories_produced_per_day: float = 0.0
o2_surplus_mol_s: float = 0.0
# Capacity Limits
human_nodes_food_limit: float = 0.0
human_nodes_o2_limit: float = 0.0
active_human_nodes: float = 0.0
system_bottleneck: str = ""
@dataclass
class NitrateChainResult:
# SMR + WGS
ch4_fed_to_smr_mol_s: float = 0.0
h2o_fed_to_smr_mol_s: float = 0.0
h2_from_smr_mol_s: float = 0.0
co2_from_smr_mol_s: float = 0.0 # recycled back to DAC loop
smr_heat_consumed_W: float = 0.0 # endothermic draw on hot side
# Haber-Bosch
n2_available_mol_s: float = 0.0
h2_available_for_haber_mol_s: float = 0.0
n2_consumed_mol_s: float = 0.0
h2_consumed_haber_mol_s: float = 0.0
nh3_produced_mol_s: float = 0.0
haber_heat_liberated_W: float = 0.0
# O2 from electrolysis byproduct
o2_from_electrolysis_mol_s: float = 0.0
o2_surplus_mol_s: float = 0.0
# Ostwald (NH3 → HNO3)
nh3_fed_ostwald_mol_s: float = 0.0
o2_consumed_ostwald_mol_s: float = 0.0
hno3_produced_mol_s: float = 0.0
h2o_produced_ostwald_mol_s: float = 0.0
ostwald_heat_liberated_W: float = 0.0
# Nitrate product
nitrate_product: str = "NH4NO3"
nitrate_mol_s: float = 0.0
nh3_used_neutralisation_mol_s: float = 0.0
# Heat balance
net_heat_to_hot_side_W: float = 0.0 # positive = adds to Stirling hot side
loop_still_closed: bool = False
notes: List[str] = field(default_factory=list) # type: ignore[assignment]
# ── simulator ────────────────────────────────────────────────────────────────
def _air_co2_mol_per_s(cfg: LoopConfig) -> float:
"""Moles of CO2 per second in an air stream based on environmental PV=nRT."""
vol_flow_m3_s = cfg.dac_air_flow_L_s * 1e-3
mol_air_per_s = (cfg.atmospheric_pressure_Pa * vol_flow_m3_s) / (R_gas * cfg.ambient_temp_K)
return mol_air_per_s * cfg.co2_mole_fraction
def _air_h2o_mol_per_s(cfg: LoopConfig) -> float:
"""Approximate moles of H2O vapour per second based on temperature and RH."""
# simple Antoine eq for water vapor pressure over water
T_C = cfg.ambient_temp_K - 273.15
if T_C > 0:
# standard Antoine for 0-100 C
P_sat = 10 ** (8.07131 - 1730.63 / (233.426 + T_C)) * 133.322 # mmHg to Pa
else:
P_sat = 0.0 # Ignore sublimation/ice partial pressure for simplicity at this scale
P_atm = cfg.atmospheric_pressure_Pa
x_h2o = (cfg.relative_humidity * P_sat) / P_atm
vol_flow_m3_s = cfg.dac_air_flow_L_s * 1e-3
mol_air_per_s = (cfg.atmospheric_pressure_Pa * vol_flow_m3_s) / (R_gas * cfg.ambient_temp_K)
return mol_air_per_s * x_h2o
def run_simulation(cfg: LoopConfig) -> LoopResult:
res = LoopResult()
notes = res.energy_balance_notes
# ── 1. solar heat input ─────────────────────────────────────────────────
solar_thermal = (
cfg.solar_irradiance_W_m2
* cfg.collector_area_m2
* cfg.collector_optical_efficiency
)
solar_pv = (
cfg.solar_irradiance_W_m2
* cfg.pv_area_m2
* cfg.pv_efficiency
)
res.solar_thermal_W = solar_thermal
res.solar_pv_W = solar_pv
cfg.stirling_heat_input_W = solar_thermal # will grow after Sabatier recovery
# ── 2. Stirling cycle ───────────────────────────────────────────────────
T_h = cfg.T_hot_K
T_c = cfg.T_cold_K
carnot = 1.0 - T_c / T_h
actual_eff = carnot * cfg.stirling_efficiency_fraction
res.carnot_efficiency = round(carnot, 4)
res.actual_efficiency = round(actual_eff, 4)
# Initial PV work before Sabatier heat recovery
pv_work_W = cfg.stirling_heat_input_W * actual_eff
heat_rejected_W = cfg.stirling_heat_input_W - pv_work_W
res.stirling_heat_rejected_W = heat_rejected_W
# ── 3. Thermopile ───────────────────────────────────────────────────────
delta_T = T_h - T_c
res.thermopile_delta_T = delta_T
V_oc = cfg.thermopile_pairs * SEEBECK_V_PER_K * delta_T
res.thermopile_open_circuit_V = round(V_oc, 2)
# Max power transfer: matched load
R_int = cfg.thermopile_internal_resistance_ohm
R_load = cfg.thermopile_load_resistance_ohm
I_tp = V_oc / (R_int + R_load)
tp_power = I_tp ** 2 * R_load
res.thermopile_W = round(tp_power, 3)
total_electrical_in = solar_pv + tp_power
notes.append(f"Electrical available (PV + thermopile): {total_electrical_in:.2f} W")
# ── 4. DAC — cold side sorbent bed ──────────────────────────────────────
co2_in_mol_s = _air_co2_mol_per_s(cfg)
co2_captured = co2_in_mol_s * cfg.dac_capture_efficiency
res.co2_captured_mol_s = co2_captured
# Sorbent regen uses cold-side waste heat (TSA cycle)
dac_regen_W = co2_captured * DH_dac_regen_kJ * 1000.0 # W = mol/s * J/mol
res.dac_regen_power_W = round(dac_regen_W, 3)
if dac_regen_W > heat_rejected_W:
notes.append(
f"WARNING: DAC regen ({dac_regen_W:.1f} W) exceeds cold-side heat "
f"({heat_rejected_W:.1f} W) — reduce flow or add heat exchanger stages"
)
else:
notes.append(f"DAC regen ({dac_regen_W:.1f} W) covered by cold-side heat rejection OK")
# ── 5. H2O collection + electrolysis ────────────────────────────────────
h2o_available = _air_h2o_mol_per_s(cfg)
h2o_collected = h2o_available * cfg.h2o_collection_efficiency
res.h2o_collected_mol_s = h2o_collected
# Electrolysis powered by total_electrical_in
# Power to electrolyze: P = n_H2 * DH_elec / faradaic_eff
# → n_H2 = P * faradaic_eff / DH_elec
dh_elec_J = DH_electrolysis_kJ * 1000.0
h2_from_elec = (total_electrical_in * cfg.electrolysis_faradaic_efficiency) / dh_elec_J
# Cap by available H2O
h2_produced = min(h2_from_elec, h2o_collected)
res.h2_produced_mol_s = h2_produced
res.electrolysis_power_W = round(
h2_produced * dh_elec_J / cfg.electrolysis_faradaic_efficiency, 3
)
res.electrolysis_source = "solar_pv + thermopile"
notes.append(
f"H2 production: {h2_produced*1e6:.2f} µmol/s "
f"(limited by {'H2O supply' if h2_produced == h2o_collected else 'electrical power'})"
)
# ── 6. Sabatier reactor — hot zone ──────────────────────────────────────
# CO2 + 4H2 → CH4 + 2H2O
# stoichiometric H2 needed for all captured CO2
h2_needed_for_co2 = co2_captured * 4.0
if h2_produced < h2_needed_for_co2:
# H2-limited: consume all H2, fraction of CO2
h2_consumed = h2_produced
co2_consumed = h2_consumed / 4.0
co2_surplus = co2_captured - co2_consumed
h2_surplus = 0.0
else:
# CO2-limited
co2_consumed = co2_captured * cfg.sabatier_conversion_efficiency
h2_consumed = co2_consumed * 4.0
h2_surplus = h2_produced - h2_consumed
co2_surplus = co2_captured * (1.0 - cfg.sabatier_conversion_efficiency)
ch4_produced = co2_consumed # 1:1 molar
sabatier_heat_W = co2_consumed * abs(DH_sabatier_kJ) * 1000.0
heat_recovery_W = sabatier_heat_W * cfg.sabatier_heat_recovery_fraction
res.ch4_produced_mol_s = ch4_produced
res.sabatier_heat_liberated_W = round(sabatier_heat_W, 3)
res.sabatier_heat_recovery_W = round(heat_recovery_W, 3)
res.h2_consumed_mol_s = h2_consumed
res.co2_consumed_mol_s = co2_consumed
res.h2_surplus_mol_s = round(h2_surplus, 9)
res.co2_surplus_mol_s = round(co2_surplus, 9)
# ── 6b. NITRATE CHAIN (SMR, Haber, Ostwald) ─────────────────────────────
# Treat the system as an unfolded 3D device tracking atomic traces (C, H, N, O)
nitrate_cfg = NitrateChainConfig()
nc = NitrateChainResult()
# ── Methane diversion
ch4_surplus = ch4_produced * nitrate_cfg.ch4_surplus_fraction
ch4_to_compressor = ch4_produced - ch4_surplus
# SMR trace: C and H nodes out of loop. CH4 + 2H2O -> CO2 + 4H2
ch4_to_smr = ch4_surplus * nitrate_cfg.ch4_smr_fraction
nc.ch4_fed_to_smr_mol_s = ch4_to_smr
nc.h2o_fed_to_smr_mol_s = ch4_to_smr * 2.0
nc.co2_from_smr_mol_s = ch4_to_smr
nc.h2_from_smr_mol_s = ch4_to_smr * 4.0
nc.smr_heat_consumed_W = ch4_to_smr * DH_smr_net_kJ * 1000.0 # Endothermic
# Nitrogen & Hydrogen traces connecting into Haber
h2_avail = h2_surplus + nc.h2_from_smr_mol_s
nc.h2_available_for_haber_mol_s = h2_avail
# N2 from free air trace
vol_flow_m3_s = cfg.dac_air_flow_L_s * 1e-3
mol_air_per_s = (cfg.atmospheric_pressure_Pa * vol_flow_m3_s) / (R_gas * cfg.ambient_temp_K)
n2_avail = mol_air_per_s * cfg.n2_mole_fraction
nc.n2_available_mol_s = n2_avail
# Haber node resolution: N2 + 3H2 -> 2NH3
n2_needed_for_h2 = h2_avail / 3.0
n2_reacted = min(n2_avail, n2_needed_for_h2) * nitrate_cfg.haber_n2_conversion
nc.h2_consumed_haber_mol_s = n2_reacted * 3.0
nc.n2_consumed_mol_s = n2_reacted
nc.nh3_produced_mol_s = n2_reacted * 2.0
nc.haber_heat_liberated_W = n2_reacted * abs(DH_haber_kJ_per_n2) * 1000.0
# Ostwald trace linking NH3 and Oxygen (from Electrolysis byproduct)
nh3_to_ostwald = nc.nh3_produced_mol_s * 0.5 # Split 50/50 for NH4NO3
nc.o2_from_electrolysis_mol_s = h2_produced * 0.5
nc.nh3_fed_ostwald_mol_s = nh3_to_ostwald
nc.hno3_produced_mol_s = nh3_to_ostwald * nitrate_cfg.ostwald_nh3_conversion
nc.ostwald_heat_liberated_W = nh3_to_ostwald * abs(DH_ostwald_kJ_per_nh3) * 1000.0
# Neutralisation node: Crossover mapping to topsoil salt
nc.nitrate_mol_s = min(nc.nh3_produced_mol_s - nh3_to_ostwald, nc.hno3_produced_mol_s)
nc.nh3_used_neutralisation_mol_s = nc.nitrate_mol_s
# Heat aggregate node linking back to Stirling loop
nc.net_heat_to_hot_side_W = (nc.haber_heat_liberated_W + nc.ostwald_heat_liberated_W) - nc.smr_heat_consumed_W
res.nitrate_chain = nc
# ── 6c. BIOLOGICAL LOAD (Consumer Packets routing) ──────────────────────
bio_cfg = BiologicalLoadConfig()
bio = BiologicalLoadResult()
# Each mol of NH4NO3 gives 2 mols of Nitrogen atoms for biomass
bio.n_routed_to_biomass_mol_s = nc.nitrate_mol_s * 2.0
# 1 mol N = 14.0067 grams.
n_kg_s = bio.n_routed_to_biomass_mol_s * 0.0140067
# If biomass is 1.5% N, total biomass = n_kg / 0.015
bio.biomass_kg_s = n_kg_s / bio_cfg.crop_n_fraction
# Calories generated per day = biomass_kg_s * (seconds in day) * kcal_per_kg
bio.calories_produced_per_day = bio.biomass_kg_s * 86400.0 * bio_cfg.biomass_kcal_per_kg
bio.human_nodes_food_limit = bio.calories_produced_per_day / bio_cfg.human_kcal_per_day
# O2 routing: electrolysis O2 - ostwald O2
bio.o2_surplus_mol_s = nc.o2_from_electrolysis_mol_s - nc.o2_consumed_ostwald_mol_s
if bio.o2_surplus_mol_s < 0:
bio.o2_surplus_mol_s = 0.0
o2_mol_day = bio.o2_surplus_mol_s * 86400.0
bio.human_nodes_o2_limit = o2_mol_day / bio_cfg.human_o2_mol_per_day
bio.active_human_nodes = min(bio.human_nodes_food_limit, bio.human_nodes_o2_limit)
bio.system_bottleneck = "Food" if bio.human_nodes_food_limit < bio.human_nodes_o2_limit else ("O2" if bio.human_nodes_o2_limit < bio.human_nodes_food_limit else "Balanced")
res.bio_load = bio
# Feed Sabatier and Nitrate node heat deltas back to Stirling hot side
total_heat_in = solar_thermal + heat_recovery_W + nc.net_heat_to_hot_side_W
res.total_heat_in_W = round(total_heat_in, 3)
pv_work_W = total_heat_in * actual_eff
stirling_electrical = pv_work_W * cfg.linear_alternator_efficiency
res.stirling_pv_work_W = round(pv_work_W, 3)
res.stirling_electrical_W = round(stirling_electrical, 3)
notes.append(
f"Thermochemical recovery (Sabatier + Nitrate chain) shifts hot side energy by "
f"{(heat_recovery_W + nc.net_heat_to_hot_side_W):.2f} W → "
f"AC flow wire resolved: {stirling_electrical:.2f} W"
)
# ── 7. Compressor → supercritical product ───────────────────────────────
if cfg.product_mode == "CH4":
mol_s = ch4_to_compressor
T_crit = T_crit_ch4
P_target = P_crit_ch4 * 1.3 # 30% above critical
else:
mol_s = co2_captured - co2_consumed + co2_surplus
T_crit = T_crit_co2
P_target = P_crit_co2 * 1.3
res.product_mol_s = mol_s
# Isothermal compression work estimate: W = n·R·T·ln(P2/P1)
P_in = 101325.0 # Pa
if mol_s > 0 and P_target > P_in:
W_ideal = mol_s * R_gas * T_std * math.log(P_target / P_in)
W_actual = W_ideal / cfg.compressor_isentropic_efficiency
else:
W_actual = 0.0
compressor_electrical = 0.0 if not W_actual else W_actual / cfg.compressor_motor_efficiency
res.compressor_shaft_input_W = round(W_actual, 3)
res.compressor_electrical_input_W = round(compressor_electrical, 3)
res.product_pressure_Pa = P_target
res.product_temp_K = T_std # post-intercooling
# Check supercritical condition
near_critical = abs(T_std - T_crit) < 50.0 # within 50 K of critical T
above_critical_P = P_target > (P_crit_ch4 if cfg.product_mode == "CH4" else P_crit_co2)
res.supercritical = above_critical_P and near_critical
# ── 8. Loop closure (Micro-Grid Electrical Bus) ─────────────────────────
# FPSCs decouple the compressor via electrical bus
elec_consumed = res.electrolysis_power_W + compressor_electrical
elec_available = total_electrical_in + res.stirling_electrical_W
net_elec = elec_available - elec_consumed
res.net_electrical_surplus_W = round(net_elec, 3)
# Loop is closed if:
# 1. Net electrical grid is positive (covers Electrolysis + Compressor)
# 2. DAC regen covered by cold-side heat
elec_ok = net_elec >= 0
dac_ok = dac_regen_W <= heat_rejected_W
res.loop_closed = elec_ok and dac_ok
if not elec_ok:
notes.append(
f"LOOP OPEN: Micro-Grid electrical deficit {abs(net_elec):.2f} W — "
f"increase collector capacity or reduce conversion load"
)
if not dac_ok:
notes.append(
f"LOOP OPEN: DAC heat deficit {dac_regen_W - heat_rejected_W:.2f} W — "
f"reduce air flow or add auxiliary heat exchanger"
)
if res.loop_closed:
notes.append(
"LOOP CLOSED: all electrical and DAC thermal budgets successfully routed "
"via Free-Piston electrical grid."
)
return res
# ── output renderers ────────────────────────────────────────────────────────
def render_report(cfg: LoopConfig, res: LoopResult) -> None:
W = 58
def row(label: str, value: str) -> str:
return f" {label:<36} {value}"
print(f"\n{'═' * W}")
print(f" STIRLING-DAC-SABATIER CLOSED LOOP {'═' * (W - 36)}")
print(" [ SYSTEM STATE: WAVEFORM COLLAPSED -> GROUND STATE ]")
print(" [ N-TRACE EXCITATION: FPSC ACOUSTIC SONICATION ACTIVE ]")
if res.nitrate_chain and res.nitrate_chain.nitrate_mol_s > 0:
nc = res.nitrate_chain
print(f"\n {'── NITRATE TRACE WIRES (Atomic topological graph)':}")
print(row("C-Trace: CH4 bypassed to SMR", f"{nc.ch4_fed_to_smr_mol_s*1e6:.4f} µmol/s"))
print(row("H-Trace: H2 harvested off SMR", f"{nc.h2_from_smr_mol_s*1e6:.4f} µmol/s"))
print(row("N-Trace: N2 entrained via air", f"{nc.n2_consumed_mol_s*1e6:.4f} µmol/s"))
print(row("O-Trace: O2 wire from PV H2O lysis", f"{nc.o2_from_electrolysis_mol_s*1e6:.4f} µmol/s"))
print(row("Node Synthesis: NH3 stabilized", f"{nc.nh3_produced_mol_s*1e6:.4f} µmol/s"))
print(row("Node Synthesis: HNO3 derived", f"{nc.hno3_produced_mol_s*1e6:.4f} µmol/s"))
print(row("Terminal: NH4NO3 soil precursor", f"{nc.nitrate_mol_s*1e6:.4f} µmol/s"))
tons_per_year = (nc.nitrate_mol_s * 80.043 * 3600 * 24 * 365.25) / 1e6
print(row("Topsoil matrix potential (1 unit/yr)", f"{tons_per_year:.4f} metric tons/year"))
print(row("Latent thermal bias via reaction", f"{nc.net_heat_to_hot_side_W:+.3f} W (to hot zone loop)"))
if res.bio_load and res.nitrate_chain and res.nitrate_chain.nitrate_mol_s > 0:
bio = res.bio_load
print(f"\n {'── BIOLOGICAL LOAD (Consumer Nodes Routing)':}")
print(row("Crop Nitrogen Matrix (N atoms)", f"{bio.n_routed_to_biomass_mol_s*1e6:.4f} µmol/s"))
print(row("Biomass Rendered (1.5% N matrix)", f"{bio.biomass_kg_s*3600*24:.3f} kg/day"))
print(row("Caloric Payload Generated", f"{bio.calories_produced_per_day:.0f} kcal/day"))
print(row("O2 Byproduct Matrix (Surplus)", f"{bio.o2_surplus_mol_s*1e6:.4f} µmol/s"))
print(row("Max Human Nodes (Food Constraint)", f"{bio.human_nodes_food_limit:.3f} souls"))
print(row("Max Human Nodes (O2 Constraint)", f"{bio.human_nodes_o2_limit:.3f} souls"))
print(row("Sustainable Human Matrix Load", f"{bio.active_human_nodes:.3f} Active Nodes (Bottleneck: {bio.system_bottleneck})"))
print(f"{'═' * W}")
print(f"\n {'── ENERGY INPUTS':}")
print(row("Solar thermal (collector)", f"{res.solar_thermal_W:.2f} W"))
print(row("Solar PV (electrolysis)", f"{res.solar_pv_W:.2f} W"))
print(row("Thermopile (Stirling ΔT)", f"{res.thermopile_W:.3f} W"))
print(row("Sabatier exotherm recovery", f"{res.sabatier_heat_recovery_W:.3f} W"))
print(row("Total hot-side heat in", f"{res.total_heat_in_W:.2f} W"))
print(f"\n {'── STIRLING (NASA GRC Free-Piston Topology)':}")
print(row("Carnot efficiency", f"{res.carnot_efficiency*100:.1f}%"))
print(row("Actual PV efficiency (×{:.0f}% Carnot)".format(
cfg.stirling_efficiency_fraction * 100), f"{res.actual_efficiency*100:.1f}%"))
print(row("PV thermodynamic work", f"{res.stirling_pv_work_W:.2f} W"))
print(row("Linear alternator output (AC)", f"{res.stirling_electrical_W:.2f} W"))
print(row("Heat rejected (radiator sink)", f"{res.stirling_heat_rejected_W:.2f} W"))
print(f"\n {'── THERMOPILE':}")
print(row("ΔT hot/cold", f"{res.thermopile_delta_T:.0f} K"))
print(row("Open-circuit voltage", f"{res.thermopile_open_circuit_V:.1f} V"))
print(row("Power (matched load)", f"{res.thermopile_W:.3f} W"))
print(f"\n {'── DIRECT AIR CAPTURE (cold side)':}")
print(row("CO2 captured", f"{res.co2_captured_mol_s*1e6:.3f} µmol/s"))
print(row("Regen heat needed", f"{res.dac_regen_power_W:.3f} W"))
print(row("Cold-side heat available", f"{res.stirling_heat_rejected_W:.2f} W"))
print(f"\n {'── H2O + ELECTROLYSIS':}")
print(row("H2O collected (condensate)", f"{res.h2o_collected_mol_s*1e6:.3f} µmol/s"))
print(row("H2 produced", f"{res.h2_produced_mol_s*1e6:.3f} µmol/s"))
print(row("Electrolysis power consumed", f"{res.electrolysis_power_W:.3f} W"))
print(row("Source", res.electrolysis_source))
print(f"\n {'── SABATIER REACTOR (hot zone)':}")
print(row("CO2 consumed", f"{res.co2_consumed_mol_s*1e6:.4f} µmol/s"))
print(row("H2 consumed", f"{res.h2_consumed_mol_s*1e6:.4f} µmol/s"))
print(row("CH4 produced", f"{res.ch4_produced_mol_s*1e6:.4f} µmol/s"))
print(row("Exotherm liberated", f"{res.sabatier_heat_liberated_W:.4f} W"))
print(row(" → fed back to hot side", f"{res.sabatier_heat_recovery_W:.4f} W"))
print(row("H2 surplus", f"{res.h2_surplus_mol_s*1e9:.3f} nmol/s"))
print(row("CO2 surplus", f"{res.co2_surplus_mol_s*1e9:.3f} nmol/s"))
print(f"\n {'── COMPRESSOR → SUPERCRITICAL':}")
print(row("Product mode", cfg.product_mode))
print(row("Product flow", f"{res.product_mol_s*1e6:.4f} µmol/s"))
print(row("Target pressure", f"{res.product_pressure_Pa/1e6:.3f} MPa"))
print(row("Is supercritical?", "YES ✓" if res.supercritical else "no (sub-Pc or ΔT >50 K)"))
print(row("Shaft work required", f"{res.compressor_shaft_input_W:.4f} W"))
print(row("Motor electric overhead", f"{res.compressor_electrical_input_W:.4f} W"))
print(f"\n {'── LOOP CLOSURE (Unified Micro-Grid)':}")
print(row("Net electrical surplus", f"{res.net_electrical_surplus_W:.3f} W"))
status = "✓ CLOSED — self-sustaining" if res.loop_closed else "✗ OPEN — see notes"
print(row("Loop status", status))
print(f"\n {'── NOTES':}")
for note in res.energy_balance_notes:
for line in (note[i:i+W-4] for i in range(0, len(note), W-4)):
print(f" {line}")
print(f"\n{'═' * W}\n")
def render_one_line(_cfg: LoopConfig, res: LoopResult, delim: str) -> None:
parts = [
f"solar_thermal={res.solar_thermal_W:.1f}W",
f"solar_pv={res.solar_pv_W:.1f}W",
f"thermopile={res.thermopile_W:.2f}W",
f"carnot={res.carnot_efficiency*100:.1f}pct",
f"actual_eff={res.actual_efficiency*100:.1f}pct",
f"ac_out={res.stirling_electrical_W:.2f}W",
f"thermopile_dT={res.thermopile_delta_T:.0f}K",
f"co2_captured={res.co2_captured_mol_s*1e6:.3f}umol_s",
f"h2_prod={res.h2_produced_mol_s*1e6:.3f}umol_s",
f"ch4_prod={res.ch4_produced_mol_s*1e6:.4f}umol_s",
f"supercritical={'yes' if res.supercritical else 'no'}",
f"elec_surplus={res.net_electrical_surplus_W:.3f}W",
]
if res.nitrate_chain:
parts.append(f"nh4no3_prod_umol={res.nitrate_chain.nitrate_mol_s*1e6:.4f}")
parts.append(f"loop={'CLOSED' if res.loop_closed else 'OPEN'}")
print(delim.join(parts))
# ── CLI ─────────────────────────────────────────────────────────────────────
def main() -> None:
ap = argparse.ArgumentParser(
description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
)
ap.add_argument("--t-hot", type=float, default=954.2, metavar="K",
help="Stirling hot side temperature K (default: 954.2 [Simulated Annealing Minimum])")
ap.add_argument("--t-cold", type=float, default=363.0, metavar="K",
help="Stirling cold side temperature K (default: 363 = 90°C)")
ap.add_argument("--stirling-eff-fraction", type=float, default=0.55, metavar="F",
help="Fraction of Carnot efficiency achieved (default: 0.55)")
ap.add_argument("--solar-w", type=float, default=800.0, metavar="W/m2",
help="Solar irradiance W/m² (default: 800)")
ap.add_argument("--area-m2", type=float, default=1.5, metavar="M2",
help="Solar collector area m² (default: 1.5)")
ap.add_argument("--pv-area-m2", type=float, default=0.20, metavar="M2",
help="PV panel area m² for electrolysis (default: 0.20)")
ap.add_argument("--air-flow-l-s", type=float, default=5.0, metavar="L/s",
help="DAC air flow rate L/s (default: 5.0)")
ap.add_argument("--tp-pairs", type=int, default=120, metavar="N",
help="Number of thermopile TE couples (default: 120)")
ap.add_argument("--product", choices=["CH4", "CO2"], default="CH4",
help="Supercritical product: CH4 or CO2 (default: CH4)")
ap.add_argument("--atm-pressure", type=float, default=101325.0, metavar="Pa",
help="Atmospheric pressure in Pa (default: 101325)")
ap.add_argument("--ambient-temp", type=float, default=298.15, metavar="K",
help="Ambient temperature in K (default: 298.15)")
ap.add_argument("--co2-ppm", type=float, default=420.0, metavar="PPM",
help="CO2 concentration in PPM (default: 420)")
ap.add_argument("--n2-frac", type=float, default=0.7809, metavar="FRAC",
help="N2 mole fraction (default: 0.7809)")
ap.add_argument("--o2-frac", type=float, default=0.2095, metavar="FRAC",
help="O2 mole fraction (default: 0.2095)")
ap.add_argument("--rh", type=float, default=0.50, metavar="FRAC",
help="Relative humidity fraction (default: 0.50)")
ap.add_argument("--json", action="store_true", help="Output full JSON report")
ap.add_argument("--one-line", action="store_true", help="Pipe-friendly one-line output")
ap.add_argument("--one-line-delim", default=" | ", metavar="DELIM",
help="Delimiter for --one-line (default: ' | ')")
args = ap.parse_args()
cfg = LoopConfig(
T_hot_K=args.t_hot,
T_cold_K=args.t_cold,
stirling_efficiency_fraction=args.stirling_eff_fraction,
solar_irradiance_W_m2=args.solar_w,
collector_area_m2=args.area_m2,
pv_area_m2=args.pv_area_m2,
thermopile_pairs=args.tp_pairs,
dac_air_flow_L_s=args.air_flow_l_s,
product_mode=args.product,
atmospheric_pressure_Pa=args.atm_pressure,
ambient_temp_K=args.ambient_temp,
co2_mole_fraction=args.co2_ppm / 1e6,
n2_mole_fraction=args.n2_frac,
o2_mole_fraction=args.o2_frac,
relative_humidity=args.rh,
)
res = run_simulation(cfg)
if args.json:
out: Dict[str, Any] = {"config": asdict(cfg), "result": asdict(res)}
print(json.dumps(out, indent=2))
elif args.one_line:
render_one_line(cfg, res, args.one_line_delim)
else:
render_report(cfg, res)
if __name__ == "__main__":
main()
Reference in a new issue View git blame Copy permalink