# Thermodynamics

This document provides comprehensive documentation for thermodynamic models, property calculations, and databases available in Difflow.

---

## Overview

Difflow provides two levels of thermodynamic modeling:

| Model | Accuracy | Speed | Best For |
|-------|----------|-------|----------|
| **Ideal** | Moderate | Fast | Preliminary design, ideal mixtures |
| **Cubic EOS** | Good | Moderate | Non-ideal gases, high pressure |

All thermodynamic calculations are fully differentiable with JAX, enabling gradient-based optimization.

---

(ideal-thermodynamics)=
## Ideal Thermodynamics

**Location**: `difflow/thermo.py`

(speciesdata)=
### SpeciesData

The fundamental data structure for storing species properties.

```python
from typing import NamedTuple

class SpeciesData(NamedTuple):
    name: str              # Species name
    MW: float              # Molecular weight (g/mol)
    Cp_coeffs: tuple       # Heat capacity polynomial (a, b, c, d)
    Hvap_coeffs: tuple     # Heat of vaporization (A, n, Tc)
    antoine_coeffs: tuple  # Antoine equation (A, B, C)
    Hf: float              # Heat of formation (J/mol)
    Tref: float            # Reference temperature (K)
```

#### Heat Capacity Polynomial

$$C_p(T) = a + bT + cT^2 + dT^3$$

Where:
- $C_p$: Heat capacity at constant pressure (J/mol/K)
- $T$: Temperature (K)
- $a, b, c, d$: Polynomial coefficients

**Units**:
- $a$: J/mol/K
- $b$: J/mol/K²
- $c$: J/mol/K³
- $d$: J/mol/K⁴

#### Antoine Equation (Vapor Pressure)

$$\log_{10}(P^{sat}) = A - \frac{B}{T + C}$$

Where:
- $P^{sat}$: Saturation pressure (Pa)
- $T$: Temperature (K)
- $A, B, C$: Antoine coefficients

**Note**: Different sources use different forms. Difflow uses:
- Pressure in Pa
- Temperature in K

#### Watson Correlation (Heat of Vaporization)

$$\Delta H_{vap}(T) = A \left(1 - \frac{T}{T_c}\right)^n$$

Where:
- $\Delta H_{vap}$: Heat of vaporization (J/mol)
- $T_c$: Critical temperature (K)
- $A, n$: Watson correlation parameters

**Example**:

```python
from difflow.thermo import SpeciesData

# Define ethanol
ethanol = SpeciesData(
    name='ethanol',
    MW=46.07,                          # g/mol
    Cp_coeffs=(9.014, 0.2141, -8.39e-5, 1.373e-8),  # J/mol/K
    Hvap_coeffs=(50430.0, 0.4475, 513.9),           # Watson params
    antoine_coeffs=(10.8095, 1592.86, -46.95),      # P in Pa, T in K
    Hf=-277690.0,                      # J/mol
    Tref=298.15                        # K
)
```

---

(idealthermo-class)=
### IdealThermo Class

The main class for ideal thermodynamic calculations.

<!-- doc-test: skip: signature listing, redefines the class -->
```python
from difflow.thermo import IdealThermo

class IdealThermo:
    def __init__(self, species_data: dict[str, SpeciesData]):
        """
        Initialize ideal thermodynamic model.

        Args:
            species_data: Dictionary mapping species names to SpeciesData
        """
```

#### Initialization

```python
from difflow.thermo import IdealThermo
from difflow.database import get_species_data

# SpeciesData(name, MW, Cp_coeffs, Hvap_coeffs, antoine_coeffs, ...) per species;
# here taken from the built-in database
species_data = {s: get_species_data(s)
                for s in ['methanol', 'ethanol', 'water', 'dimethyl_ether']}

thermo = IdealThermo(species_data)
```

---

(property-calculations)=
### Property Calculations

#### Heat Capacity

```python
# Single species
Cp = thermo.Cp('ethanol', T=350.0)  # J/mol/K

# Mixture (molar average)
Cp_mix = thermo.Cp_mix(mole_fracs={'ethanol': 0.4, 'water': 0.6}, T=350.0)
```

**Equations**:

$$C_p^{pure}(T) = a + bT + cT^2 + dT^3$$

$$C_p^{mix} = \sum_i x_i C_{p,i}$$

#### Enthalpy

```python
# Pure species enthalpy relative to reference
H = thermo.H_pure('ethanol', T=400.0, phase='liquid')  # J/mol

# Stream enthalpy
H_flow = thermo.stream_enthalpy(
    flows={'ethanol': 1.0, 'water': 2.0},  # mol/s
    T=350.0,
    phase='liquid'
)  # W (J/s)
```

**Equations**:

$$H(T) = H_f + \int_{T_{ref}}^T C_p \, dT$$

$$H(T) = H_f + a(T - T_{ref}) + \frac{b}{2}(T^2 - T_{ref}^2) + \frac{c}{3}(T^3 - T_{ref}^3) + \frac{d}{4}(T^4 - T_{ref}^4)$$

For vapor phase, add heat of vaporization:

$$H^V(T) = H^L(T) + \Delta H_{vap}(T)$$

#### Saturation Pressure

```python
P_sat = thermo.Psat('ethanol', T=350.0)  # Pa
```

**Equation** (Antoine):

$$P^{sat} = 10^{A - B/(T + C)}$$

#### Heat of Vaporization

```python
Hvap = thermo.Hvap('ethanol', T=350.0)  # J/mol
```

**Equation** (Watson correlation):

$$\Delta H_{vap} = A \left(1 - \frac{T}{T_c}\right)^n$$

#### K-Values (Vapor-Liquid Equilibrium)

```python
# Single species K-value (Raoult's law)
K = thermo.K_value('ethanol', T=350.0, P=101325.0)

# All species K-values
K_values = thermo.K_values(T=350.0, P=101325.0)  # dict
```

**Equation** (Raoult's Law):

$$K_i = \frac{y_i}{x_i} = \frac{P_i^{sat}(T)}{P}$$

**Assumptions**:
- Ideal liquid mixture (activity coefficient = 1)
- Ideal gas phase (fugacity coefficient = 1)
- Valid for low pressures and similar molecules

`K_values` and `K_values_array` also accept the liquid and vapor compositions
`x` and `y` and ignore them — Raoult K-values are a function of $(T, P)$ alone.
They are in the signature so that a unit written against this interface (the
rigorous [`DistillationColumn`](unit-operations-chemical.md), for one) can pass
its stage compositions unconditionally and run unchanged on a
[`CubicThermo`](#cubic-thermo-k-values), whose K-values do depend on them.
`IdealThermo.K_depends_on_composition` is `False` and `CubicThermo`'s is `True`,
for callers that want to skip a composition iteration they do not need.

#### Bubble Point and Dew Point

```python
# Bubble pressure at given T and liquid composition
P_bubble = thermo.bubble_pressure(x={'ethanol': 0.4, 'water': 0.6}, T=350.0)

# Dew pressure at given T and vapor composition
P_dew = thermo.dew_pressure(y={'ethanol': 0.4, 'water': 0.6}, T=350.0)
```

**Equations**:

Bubble point: $P = \sum_i x_i P_i^{sat}$

Dew point: $\frac{1}{P} = \sum_i \frac{y_i}{P_i^{sat}}$

---

(cubic-equations-of-state)=
## Cubic Equations of State

**Location**: `difflow/eos.py`

Cubic equations of state provide more accurate thermodynamic predictions for non-ideal systems, especially at high pressures.

### Critical Properties

<!-- doc-test: skip: signature listing, redefines the class -->
```python
from difflow.eos import CriticalProperties

class CriticalProperties(NamedTuple):
    name: str          # Species name
    Tc: float          # Critical temperature (K)
    Pc: float          # Critical pressure (Pa)
    omega: float       # Acentric factor
    MW: float          # Molecular weight (g/mol)
```

**Acentric Factor** ($\omega$):

$$\omega = -\log_{10}\left(\frac{P^{sat}(T_r=0.7)}{P_c}\right) - 1$$

Measures deviation from simple fluid behavior:
- $\omega \approx 0$: Spherical molecules (Ar, Kr)
- $\omega > 0$: Non-spherical or polar molecules

(peng-robinson-eos)=
### Peng-Robinson EOS

The Peng-Robinson equation of state (1976) is widely used for hydrocarbon systems.

```python
from difflow.eos import PengRobinson, CriticalProperties

# Initialize with species critical properties
critical_props = {
    'methane': CriticalProperties('methane', 190.6, 4.6e6, 0.011, 16.04),
    'ethane': CriticalProperties('ethane', 305.4, 4.88e6, 0.099, 30.07),
}

pr = PengRobinson(critical_props)
```

#### Equations

**Equation of State**:

$$P = \frac{RT}{V - b} - \frac{a(T)}{V^2 + 2bV - b^2}$$

Or in terms of compressibility factor $Z = PV/RT$:

$$Z^3 - (1-B)Z^2 + (A - 3B^2 - 2B)Z - (AB - B^2 - B^3) = 0$$

Where:
- $A = aP/(R^2T^2)$
- $B = bP/(RT)$

**Parameters**:

$$a(T) = a_c \cdot \alpha(T)$$

$$a_c = 0.45724 \frac{R^2 T_c^2}{P_c}$$

$$b = 0.07780 \frac{RT_c}{P_c}$$

$$\alpha(T) = \left[1 + \kappa\left(1 - \sqrt{T/T_c}\right)\right]^2$$

$$\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2$$

**Mixing Rules** (van der Waals one-fluid):

$$a_{mix} = \sum_i \sum_j y_i y_j \sqrt{a_i a_j}(1 - k_{ij})$$

$$b_{mix} = \sum_i y_i b_i$$

Where $k_{ij}$ is the binary interaction parameter (default = 0).

#### Methods

```python
import jax.numpy as jnp
from difflow.eos import flash_TP_eos

# Compositions are arrays in pr.species_order (here methane, ethane)
y = jnp.array([0.7, 0.3])

# Compressibility factor (largest root = vapor, smallest = liquid)
Z = pr.solve_Z(T=300.0, P=1e6, y=y, phase='vapor')

# Fugacity coefficients of every species
phi = pr.fugacity_coefficient(T=300.0, P=1e6, y=y, phase='vapor')

# K-values from fugacity, given liquid and vapor compositions
x = jnp.array([0.4, 0.6])
K = pr.K_values(T=250.0, P=2e6, x=x, y=y)

# VLE flash calculation: vapor fraction, liquid and vapor compositions
V_frac, x, y = flash_TP_eos(pr, z=jnp.array([0.5, 0.5]), T=250.0, P=2e6)
```

#### Fugacity Calculation

$$\ln \phi_i = \frac{b_i}{b_{mix}}(Z - 1) - \ln(Z - B) - \frac{A}{2\sqrt{2}B}\left(\frac{2\sum_j y_j a_{ij}}{a_{mix}} - \frac{b_i}{b_{mix}}\right)\ln\left(\frac{Z + (1+\sqrt{2})B}{Z + (1-\sqrt{2})B}\right)$$

**Equilibrium Condition**:

$$f_i^V = f_i^L$$

$$y_i \phi_i^V P = x_i \phi_i^L P$$

$$K_i = \frac{y_i}{x_i} = \frac{\phi_i^L}{\phi_i^V}$$

---

(soave-redlich-kwong-eos)=
### Soave-Redlich-Kwong EOS

The SRK equation (1972) is another popular cubic EOS.

```python
from difflow.eos import SRK

srk = SRK(critical_props)
```

#### Equations

**Equation of State**:

$$P = \frac{RT}{V - b} - \frac{a(T)}{V(V + b)}$$

**Parameters**:

$$a_c = 0.42748 \frac{R^2 T_c^2}{P_c}$$

$$b = 0.08664 \frac{RT_c}{P_c}$$

$$\alpha(T) = \left[1 + m\left(1 - \sqrt{T/T_c}\right)\right]^2$$

$$m = 0.480 + 1.574\omega - 0.176\omega^2$$

### Comparison: PR vs SRK

| Property | Peng-Robinson | SRK |
|----------|---------------|-----|
| Liquid density | Better | Less accurate |
| Vapor pressure | Good | Good |
| Near critical | Better | Good |
| Polar compounds | Limited | Limited |
| Parameters | $\Omega_a = 0.45724$ | $\Omega_a = 0.42748$ |
| | $\Omega_b = 0.07780$ | $\Omega_b = 0.08664$ |

---

(flash-calculations)=
### Flash Calculations

Flash calculations determine phase split at specified T and P.

```python
import jax.numpy as jnp
from difflow.eos import PengRobinson, flash_TP_eos
from difflow.database import get_critical_props

names = ['methane', 'ethane', 'propane']
pr3 = PengRobinson({n: get_critical_props(n) for n in names})

# TP Flash using PR EOS (compositions are arrays in species order)
V_frac, x, y = flash_TP_eos(pr3, z=jnp.array([0.3, 0.3, 0.4]), T=250.0, P=1.5e6)

print(f"Vapor fraction: {float(V_frac):.3f}")
print(f"Liquid composition: {dict(zip(names, map(float, x)))}")
print(f"Vapor composition: {dict(zip(names, map(float, y)))}")
```

#### Algorithm

1. **Initial K-values** (Wilson correlation):
   $$K_i = \frac{P_{c,i}}{P} \exp\left[5.373(1 + \omega_i)(1 - T_{c,i}/T)\right]$$

2. **Rachford-Rice equation** (solve for V):
   $$f(V) = \sum_i \frac{z_i(K_i - 1)}{1 + V(K_i - 1)} = 0$$

3. **Phase compositions**:
   $$x_i = \frac{z_i}{1 + V(K_i - 1)}$$
   $$y_i = K_i x_i$$

4. **Update K-values** from fugacity:
   $$K_i^{new} = \frac{\phi_i^L}{\phi_i^V}$$

5. **Iterate** until convergence

---

(cubic-thermo-k-values)=
### CubicThermo K-Values

`CubicThermo` exposes the same K-value interface as `IdealThermo`, built from
the EOS fugacity coefficients rather than from Raoult's law:

$$K_i = \frac{\hat\phi_i^L(T, P, x)}{\hat\phi_i^V(T, P, y)}$$

```python
from difflow.thermo import CubicThermo

names = ['propane', 'n_butane']
sp = {n: get_species_data(n) for n in names}
crit = {n: get_critical_props(n) for n in names}
x = jnp.array([0.5, 0.5])
y = jnp.array([0.8, 0.2])

thermo = CubicThermo(IdealThermo(sp), PengRobinson(crit))

K = thermo.K_values_array(T=380.0, P=10e5, x=x, y=y)  # both compositions known
K = thermo.K_values_array(T=380.0, P=10e5, x=x)       # bubble-point K at x
K = thermo.K_values_array(T=380.0, P=10e5)            # no composition: Raoult
K = thermo.K_values(T=380.0, P=10e5, x=x)             # same, as a dict
```

That is what lets a unit written against `IdealThermo`'s interface — the
rigorous [`DistillationColumn`](unit-operations-chemical.md) — run on
Peng-Robinson without changing the unit.

Two properties of these K-values shape how they are used:

**They depend on composition.** $K_i$ is a fixed point, not a formula. Pass
whichever compositions you have; whichever you omit is filled in by a bounded
successive substitution ($y = Kx$ renormalised, or $x = y/K$) from a
composition-free starting estimate, which is the standard bubble- or dew-point
K calculation. Pass both when you already have a consistent pair — that is a
single evaluation with no inner loop.

**They only exist inside the two-root window.** The EOS gives a two-phase K
only where its cubic has two distinct roots at this $(T, P, x)$. Away from the
bubble point — a subcooled liquid, a superheated vapor — there is one root,
both phases take it, and $K_i$ comes back identically 1. That is the EOS
correctly reporting a single phase, but it makes $\sum_i K_i x_i - 1$ a flat
zero, which a root finder reads as converged wherever it is standing. A
bubble-point solve on these K-values therefore needs to start inside the window
and be able to retreat if a step leaves it; see the two-pass solve in
`difflow.units.distillation._bubble_T`.

With no composition at all, `CubicThermo` returns the wrapped `IdealThermo`'s
Raoult K-values rather than the EOS's own Wilson estimate. Both are
composition-free, but Antoine coefficients are fitted vapor-pressure data,
while Wilson is a two-constant fit off the critical point that for a heavy
hydrocarbon can be a hundred degrees out — far enough to start the EOS
iteration outside the window.

---

(species-database)=
## Species Database

**Location**: `difflow/database.py`

(available-species)=
### Available Species

The database contains ~100+ species with complete thermodynamic data:

#### Light Gases
- Hydrogen (H2), Helium (He), Nitrogen (N2), Oxygen (O2)
- Carbon monoxide (CO), Carbon dioxide (CO2)
- Hydrogen sulfide (H2S), Ammonia (NH3), Sulfur dioxide (SO2)

#### Alkanes (C1-C10)
- Methane, Ethane, Propane, n-Butane, i-Butane
- n-Pentane, i-Pentane, Neopentane
- n-Hexane, n-Heptane, n-Octane, n-Nonane, n-Decane

#### Alkenes
- Ethylene, Propylene
- 1-Butene, cis-2-Butene, trans-2-Butene, Isobutylene

#### Aromatics (BTEX)
- Benzene, Toluene
- o-Xylene, m-Xylene, p-Xylene
- Ethylbenzene, Styrene

#### Alcohols
- Methanol, Ethanol
- 1-Propanol, 2-Propanol (Isopropanol)
- 1-Butanol, 2-Butanol

#### Ketones and Aldehydes
- Acetone, Methyl ethyl ketone (MEK)
- Formaldehyde, Acetaldehyde

#### Carboxylic Acids
- Formic acid, Acetic acid

#### Esters
- Methyl acetate, Ethyl acetate

#### Ethers
- Diethyl ether, Dimethyl ether (DME)

#### Water
- Water (H2O)

(database-functions)=
### Database Functions

```python
from difflow.database import (
    get_species_data,
    get_critical_props,
    get_species_info,
    list_species,
    get_alkanes,
    get_btex,
    get_common_solvents
)

# Get SpeciesData for ideal thermo
ethanol_data = get_species_data('ethanol')

# Get CriticalProperties for EOS
ethanol_crit = get_critical_props('ethanol')

# Get all available information
info = get_species_info('ethanol')
print(info)

# List all available species
species_list = list_species()

# Get groups of species
alkanes = get_alkanes()  # {'methane': ..., 'ethane': ..., ...}
btex = get_btex()        # CriticalProperties for the BTEX aromatics
solvents = get_common_solvents()
```

### Database Contents Example

```python
# Methanol data in database
{
    'name': 'methanol',
    'MW': 32.04,
    'Tc': 512.6,  # K
    'Pc': 8.09e6,  # Pa
    'omega': 0.566,
    'Cp_coeffs': (21.15, 7.092e-2, 2.587e-5, -2.852e-8),
    'Hvap_coeffs': (45050.0, 0.4065, 512.6),
    'antoine_coeffs': (10.2044, 1582.91, -33.50),
    'Hf': -200940.0,  # J/mol
    'Tref': 298.15
}
```

---

(cantera-import)=
## Cantera Import

**Location**: `difflow/cantera_import.py`

Import thermodynamic data from Cantera YAML mechanism files without requiring Cantera installation.

(importing-mechanisms)=
### Importing Mechanisms

<!-- doc-test: skip: needs a Cantera YAML mechanism file (gri30.yaml) -->
```python
from difflow.cantera_import import (
    import_species_data,
    import_critical_props,
    import_reactions,
    load_mechanism,
    list_available_species
)

# List available species in a Cantera file
species = list_available_species('gri30.yaml')

# Import species data for ideal thermo
species_data = import_species_data(
    'gri30.yaml',
    species_list=['CH4', 'O2', 'CO2', 'H2O']
)

# Import reactions with Arrhenius kinetics
reactions = import_reactions('gri30.yaml')

# Load complete mechanism
mechanism = load_mechanism('gri30.yaml')
```

(data-conversion)=
### Data Conversion

Cantera uses NASA polynomial format for thermodynamic properties:

#### NASA 7-Coefficient Polynomial

$$\frac{C_p}{R} = a_1 + a_2 T + a_3 T^2 + a_4 T^3 + a_5 T^4$$

$$\frac{H}{RT} = a_1 + \frac{a_2}{2}T + \frac{a_3}{3}T^2 + \frac{a_4}{4}T^3 + \frac{a_5}{5}T^4 + \frac{a_6}{T}$$

$$\frac{S}{R} = a_1 \ln T + a_2 T + \frac{a_3}{2}T^2 + \frac{a_4}{3}T^3 + \frac{a_5}{4}T^4 + a_7$$

The import function converts NASA coefficients to the simpler polynomial form used in Difflow.

### Supported Cantera Data

| Data Type | Support |
|-----------|---------|
| Thermo (NASA 7) | Full |
| Thermo (NASA 9) | Full |
| Transport | Partial |
| Reactions (Arrhenius) | Full |
| Reactions (falloff) | Partial |

---

(pyglenn-import)=
## NASA Glenn (pyglenn) Import

**Location**: `difflow/pyglenn_import.py`

Import ideal-gas thermodynamic data from the NASA Glenn (CEA) thermodynamic
database, as exposed by the [`pyglenn`](https://github.com/ProfLeao/pyglenn)
package (~2030 species, NASA-9 polynomials). `pyglenn` is an **optional**
dependency:

```bash
pip install pyglenn          # or:  pip install "difflow[pyglenn]"
```

Unlike the Cantera importer (which parses a YAML file), this adapter talks to
pyglenn's `ThermochemicalCalculator` at runtime. Because its
`import_species_data` / `list_available_species` names mirror the Cantera ones,
it is exposed as a **namespace** rather than flattened into `difflow`:

<!-- doc-test: skip: needs the optional pyglenn package -->
```python
from difflow.pyglenn_import import import_species_data, list_available_species
from difflow.thermo import IdealThermo

# Find species records (id, name, phase, molecular_weight, ...)
list_available_species("CO2")

# Import ideal-gas SpeciesData for a set of species
species_data = import_species_data(["O2", "CO2", "H2O"])
thermo = IdealThermo(species_data)
```

### What is (and is not) imported

| difflow `SpeciesData` field | Source in pyglenn |
|-----------------------------|-------------------|
| `Cp_coeffs` | Cubic **least-squares fit** of pyglenn's `Cp(T)` over `T_fit_range` (default 300–1000 K) |
| `MW` | `molecular_weight` |
| `Hf` | `heat_of_formation_298K` |
| `Hvap_coeffs`, `antoine_coeffs` | **Not in NASA Glenn data** — filled with neutral placeholders, or estimated from an optional `boiling_points` / `critical_temps` you pass |

The NASA-9 form carries $1/T^2$ and $1/T$ terms that difflow's cubic
$C_p = a + bT + cT^2 + dT^3$ cannot represent exactly, so `Cp_coeffs` come from
a fit over the window you care about — set `T_fit_range` to your operating
range. Samples outside a species' valid interval (where pyglenn raises) are
dropped automatically.

```{note}
pyglenn supplies **no critical properties** (Tc, Pc, ω), so there is no
`import_critical_props` here. To build a `PengRobinson`/`SRK` or a
`CubicThermo`, pair the ideal-gas `SpeciesData` from pyglenn with
`CriticalProperties` from `difflow.database` or `difflow.cantera_import`:

    from difflow.eos import PengRobinson
    from difflow.thermo import IdealThermo, CubicThermo
    from difflow.cantera_import import import_critical_props

    sp   = import_species_data(["CH4", "CO2", "H2O"])            # ideal-gas Cp (pyglenn)
    crit = import_critical_props("gri30.yaml", ["CH4", "CO2", "H2O"])  # Tc/Pc/ω (Cantera)
    thermo = CubicThermo(IdealThermo(sp), PengRobinson(crit))    # real-gas enthalpy
```

### Supported NASA Glenn Data

| Data Type | Support |
|-----------|---------|
| Ideal-gas Cp (NASA 9) | Full (cubic fit) |
| Molecular weight | Full |
| Enthalpy of formation (298 K) | Full |
| Critical properties | Not provided by pyglenn |
| Liquid Hvap / vapor pressure | Placeholder/estimated only |

---

(dwsim-import)=
## DWSIM Import

**Location**: `difflow/dwsim_import.py`

Import compound constants and ideal-gas heat capacities from
[DWSIM](https://dwsim.org)'s thermodynamics library. DWSIM is a .NET
application, so it is reached from Python through
[`pythonnet`](https://github.com/pythonnet/pythonnet) against
`DWSIM.Thermodynamics.dll`, whose `CalculatorInterface.Calculator` is the
"DTL" calculator (an older `DWSIM.Thermodynamics.StandaloneLibrary.dll` is
used if that is what the folder holds). Checked against DWSIM 9.0.5 on Linux:

```bash
scripts/install_dwsim.sh        # DWSIM 9.0.5 .deb unpacked, .NET 8 runtime, pythonnet
export DWSIM_PATH=.../usr/local/lib/dwsim
```

DWSIM 9 is a .NET 8 build: pythonnet has to load CoreCLR (`pythonnet.load
("coreclr")`, which the importer does), not Mono, and only one CLR can be
loaded in a process, so do not `import clr` before the first import.

```{important}
Calls into DWSIM return concrete numbers through the CLR and are **not
differentiable** — JAX cannot trace through them. So, exactly like the Cantera
and pyglenn importers, this adapter uses DWSIM only as a **one-time data
source**: it reads each compound's constants and samples its ideal-gas Cp(T),
then builds difflow's own JAX-native `SpeciesData` / `CriticalProperties`.
difflow stays differentiable end to end; DWSIM is never in the gradient path.
```

Because DWSIM has critical constants (unlike pyglenn), it feeds **both**
`SpeciesData` and `CriticalProperties`, so it can build a full EOS/`CubicThermo`
on its own:

<!-- doc-test: skip: needs a DWSIM installation (.NET runtime) -->
```python
from difflow.dwsim_import import import_species_data, import_critical_props
from difflow.thermo import IdealThermo, CubicThermo
from difflow.eos import PengRobinson

from difflow.dwsim_import import DWSIMBackend, dwsim_name

names = [dwsim_name(n) for n in ("methane", "co2", "water")]   # DWSIM's names
be   = DWSIMBackend()                     # DWSIM_PATH, or dwsim_path=...; ~2 s to start
sp   = import_species_data(names, backend=be)
crit = import_critical_props(names, backend=be)
thermo = CubicThermo(IdealThermo(sp), PengRobinson(crit))
```

`DWSIM_NAMES` maps difflow's database names to DWSIM's for the refinery's
50 species, every one checked to exist in DWSIM 9.0.5 (all in its ChemSep
database). `tests/test_dwsim_import.py` (release) imports them from the real
DWSIM in a subprocess and compares Tc, Pc, omega and MW with
`difflow.database`; the differences are listed under
[Validation against DWSIM](unit-operations-refinery.md#refinery-dwsim-validation).

### What is imported (and DWSIM units)

| difflow field | DWSIM source (`ICompoundConstantProperties`) | Unit conversion |
|---------------|-----------------------------------------------|-----------------|
| `MW` | `Molar_Weight` | kg/kmol ≡ g/mol |
| `CriticalProperties.Tc/Pc/omega` | `Critical_Temperature`, `Critical_Pressure`, `Acentric_Factor` | K, Pa, — |
| `Hf` | `IG_Enthalpy_of_Formation_25C` | kJ/kg × MW → J/mol |
| `Cp_coeffs` | ideal-gas Cp via the calculator's `GetCompoundTDepProp(name, "idealGasHeatCapacity", T)` | J/mol/K, then cubic fit |
| `Hvap_coeffs`, `antoine_coeffs` | estimated from Tb/Tc | — |

```{note}
DWSIM cannot run in difflow's per-commit CI (no .NET runtime). All DWSIM
contact is isolated in `DWSIMBackend`; the import logic is backend-agnostic
and unit-tested against a fake backend every commit, and against DWSIM 9.0.5
itself in the release tier (skipped where DWSIM is not installed). For
another DWSIM build, replace `DWSIMBackend` and pass it via `backend=`.
```

---

## Usage Examples

### Complete VLE Flash with Ideal Thermo

```python
from difflow.thermo import IdealThermo
from difflow.database import get_species_data
from difflow.units.flash import Flash, FlashParams
from difflow.streams import make_stream

# Build thermo model from database
species_names = ['benzene', 'toluene', 'ethylbenzene']
species_data = {name: get_species_data(name) for name in species_names}
thermo = IdealThermo(species_data)

# Create feed stream
feed = make_stream(
    {'benzene': 0.4, 'toluene': 0.35, 'ethylbenzene': 0.25},
    T=380.0,  # K
    P=101325.0  # Pa
)

# Flash calculation
flash = Flash(FlashParams(species_order=species_names), thermo)
liquid, vapor, info = flash(feed)

print(f"Vapor fraction: {float(info['V_frac']):.3f}")
print(f"K-values: {info['K']}")
```

### High-Pressure Flash with PR EOS

```python
from difflow.eos import PengRobinson
from difflow.database import get_critical_props

# Build PR model from database
species_names = ['methane', 'ethane', 'propane', 'n_butane']
critical_props = {name: get_critical_props(name) for name in species_names}
pr = PengRobinson(critical_props)

# High-pressure flash
import jax.numpy as jnp
from difflow.eos import flash_TP_eos

z = jnp.array([0.5, 0.3, 0.15, 0.05])  # in species_names order
V_frac, x, y = flash_TP_eos(pr, z, T=250.0, P=3.0e6)  # 30 bar

print(f"Vapor fraction: {float(V_frac):.3f}")
print(f"Liquid methane: {float(x[0]):.4f}")
print(f"Vapor methane: {float(y[0]):.4f}")
```

### Sensitivity Analysis with Automatic Differentiation

```python
import jax
import jax.numpy as jnp
from difflow.thermo import IdealThermo
from difflow.database import get_species_data

# Setup
species_data = {name: get_species_data(name) for name in ['ethanol', 'water']}
thermo = IdealThermo(species_data)

# Function to differentiate
def vapor_pressure_ratio(T):
    P_eth = thermo.Psat('ethanol', T)
    P_wat = thermo.Psat('water', T)
    return P_eth / P_wat

# Gradient of vapor pressure ratio w.r.t. temperature
grad_fn = jax.grad(vapor_pressure_ratio)
sensitivity = grad_fn(350.0)
print(f"d(P_eth/P_wat)/dT at 350K: {sensitivity:.6f}")
```

---

## Best Practices

### Model Selection

| Scenario | Recommended Model |
|----------|-------------------|
| Low pressure, ideal mixtures | Ideal Thermo |
| High pressure (> 10 bar) | PR or SRK EOS |
| Hydrocarbons | PR EOS |
| Polar/non-polar mixtures | PR + Binary k_ij |
| Highly polar (water, alcohols) | Activity coefficient models* |

*Activity coefficient models (NRTL, UNIQUAC) available in LLE module.

### Temperature Ranges

- **Cp polynomial**: Valid within fitted range (typically 200-1500 K)
- **Antoine equation**: Limited range (~0.01-2 bar typically)
- **Watson correlation**: Valid T < Tc
- **Cubic EOS**: Valid for all T, better away from critical

### Numerical Stability

```python
# Use jnp.where for safe operations
def safe_K_value(P_sat, P):
    return jnp.where(P > 0, P_sat / P, 0.0)

# Avoid division by zero in LMTD
def safe_lmtd(dT1, dT2):
    ratio = dT1 / jnp.maximum(dT2, 1e-10)
    return jnp.where(
        jnp.abs(dT1 - dT2) < 1e-6,
        0.5 * (dT1 + dT2),  # Limit when dT1 ≈ dT2
        (dT1 - dT2) / jnp.log(ratio)
    )
```