Transport properties

torch_flash.transport

Transport-property correlations.

HeavyOilCSPCalibrationResult dataclass

HeavyOilCSPCalibrationResult(
    third_csp: Tensor,
    fourth_csp: Tensor,
    prediction: Tensor,
    fit: FitResult,
    sensitivity_matrix: Tensor,
    sensitivity_singular_values: Tensor,
    sensitivity_rank: int,
    sensitivity_condition_number: float,
)

Joint bounded calibration of the two Lindeloff CSP factors.

Attributes:

Name	Type	Description
third_csp, fourth_csp		Selected positive dimensionless factors in Pedersen et al. (2024), Eqs. 10.34-10.35.
prediction	Tensor	Fitted dynamic viscosities in Pa s at every calibration state.
fit	FitResult	PyTorch optimizer history and stopping diagnostics.
sensitivity_matrix	Tensor	Jacobian of log viscosity with respect to the two physical factors, with shape (states, 2).
sensitivity_singular_values	Tensor	Singular values of the sensitivity matrix.
sensitivity_rank	int	Dtype-aware numerical rank.
sensitivity_condition_number	float	Largest-to-smallest singular-value ratio, or infinity for a rank-deficient calibration.

Notes

Agreement with the supplied observations is calibration evidence, not independent validation.

HeavyOilCSPProfile dataclass

HeavyOilCSPProfile(
    viscosity: Tensor,
    liquid_composition: Tensor,
    bubble_pressure: Tensor,
    vapor_fraction: Tensor,
    converged: Tensor,
    residual_norm: Tensor,
)

Phase-aware heavy-oil viscosity states and numerical diagnostics.

Attributes:

Name	Type	Description
viscosity	Tensor	Dynamic viscosity in Pa s, one value per requested state.
liquid_composition	Tensor	Feed composition above the bubble point and equilibrium liquid composition below it, with shape (states, components).
bubble_pressure	Tensor	Bubble pressure in Pa at the corresponding state temperature.
vapor_fraction	Tensor	Equilibrium vapor fraction below the bubble point and zero above it.
converged	Tensor	Per-state phase-boundary and sub-bubble flash convergence flag.
residual_norm	Tensor	Maximum dimensionless log-fugacity residual for each sub-bubble flash; zero for homogeneous states above the bubble point.

hayduk_minhas_n_paraffin_diffusion_coefficient

hayduk_minhas_n_paraffin_diffusion_coefficient(
    temperature: Tensor,
    dynamic_viscosity: Tensor,
    solute_molar_volume: Tensor,
) -> Tensor

Estimate infinite-dilution n-paraffin diffusion in an n-paraffin liquid.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Liquid temperature in K.	required
dynamic_viscosity	Tensor	Solvent or bulk-phase dynamic viscosity in Pa s.	required
solute_molar_volume	Tensor	Diffusing solute molar volume at its normal boiling point in m3/mol.	required

Returns:

Type	Description
Tensor	Diffusion coefficient in m2/s with the broadcast input shape.

Raises:

Type	Description
InvalidStateError	If any input is nonfinite or nonpositive.

Notes

Implements Hayduk and Minhas, "Correlations for prediction of molecular diffusivities in liquids," Canadian Journal of Chemical Engineering 60 (1982), 295-299, doi:10.1002/cjce.5450600213, normal-paraffin correlation, as reproduced in Pedersen et al. (2024), Eq. 10.100. The defining equation uses solvent viscosity in cP and solute normal-boiling molar volume in cm3/mol; this SI API performs both conversions. It applies to infinite-dilution n-paraffin solutes in n-paraffin solvents, not concentrated multicomponent diffusion.

Source code in src/torch_flash/transport/diffusion.py

def hayduk_minhas_n_paraffin_diffusion_coefficient(
    temperature: Tensor,
    dynamic_viscosity: Tensor,
    solute_molar_volume: Tensor,
) -> Tensor:
    """Estimate infinite-dilution n-paraffin diffusion in an n-paraffin liquid.

    Parameters
    ----------
    temperature
        Liquid temperature in K.
    dynamic_viscosity
        Solvent or bulk-phase dynamic viscosity in Pa s.
    solute_molar_volume
        Diffusing solute molar volume at its normal boiling point in m3/mol.

    Returns
    -------
    Tensor
        Diffusion coefficient in m2/s with the broadcast input shape.

    Raises
    ------
    InvalidStateError
        If any input is nonfinite or nonpositive.

    Notes
    -----
    Implements Hayduk and Minhas, "Correlations for prediction of molecular
    diffusivities in liquids," *Canadian Journal of Chemical Engineering* 60
    (1982), 295-299, doi:10.1002/cjce.5450600213, normal-paraffin
    correlation, as reproduced in Pedersen et al. (2024), Eq. 10.100. The
    defining equation uses solvent viscosity in cP and solute normal-boiling
    molar volume in cm3/mol; this SI API performs both conversions. It applies
    to infinite-dilution n-paraffin solutes in n-paraffin solvents, not
    concentrated multicomponent diffusion.
    """
    temperature, viscosity, volume = torch.broadcast_tensors(
        temperature,
        dynamic_viscosity,
        solute_molar_volume,
    )
    if bool(
        (
            (~torch.isfinite(temperature))
            | (~torch.isfinite(viscosity))
            | (~torch.isfinite(volume))
            | (temperature <= 0.0)
            | (viscosity <= 0.0)
            | (volume <= 0.0)
        ).any()
    ):
        raise InvalidStateError("diffusion correlation requires positive finite inputs")
    viscosity_cp = 1000.0 * viscosity
    volume_cm3_mol = 1.0e6 * volume
    viscosity_exponent = 10.2 / volume_cm3_mol - 0.791
    result: Tensor = (
        13.3e-12
        * temperature.pow(1.47)
        * viscosity_cp.pow(viscosity_exponent)
        / volume_cm3_mol.pow(0.71)
    )
    return result

evaluate_heavy_oil_corresponding_states_profile

evaluate_heavy_oil_corresponding_states_profile(
    model: StateModel,
    temperature: Tensor,
    pressure: Tensor,
    feed_composition: Tensor,
    components: ComponentSet,
    *,
    third_csp: Tensor | float = 1.0,
    fourth_csp: Tensor | float = 1.0,
    initial_bubble_pressure: Tensor | float | None = None,
    tolerance: float = 1e-08,
    saturation_iterations: int = 50,
    substitution_iterations: int = 20,
    newton_iterations: int = 20,
    raise_on_failure: bool = True,
) -> HeavyOilCSPProfile

Evaluate a phase-aware live-heavy-oil viscosity pressure profile.

A liquid-vapor bubble point is solved once per unique temperature. States at or above that boundary retain the feed composition. All sub-bubble states at the same temperature are flashed simultaneously with :func:torch_flash.batched_two_phase_flash, and their equilibrium liquid compositions are passed in one batch to the heavy-oil CSP correlation.

Parameters:

Name	Type	Description	Default
model	StateModel	Homogeneous-state fugacity model used for bubble points and flashes.	required
temperature	Tensor	One-dimensional state batches in K and Pa with matching shapes.	required
pressure	Tensor	One-dimensional state batches in K and Pa with matching shapes.	required
feed_composition	Tensor	One strictly positive feed mole-fraction vector.	required
components	ComponentSet	Transport critical properties and molar masses in the same order as the model and feed.	required
third_csp	Tensor | float	Dimensionless Lindeloff calibration factors. Tensor inputs preserve gradients through the viscosity calculation.	1.0
fourth_csp	Tensor | float	Dimensionless Lindeloff calibration factors. Tensor inputs preserve gradients through the viscosity calculation.	1.0
initial_bubble_pressure	Tensor | float | None	Optional positive scalar pressure estimate in Pa reused at each temperature. Wilson initialization is used when omitted.	None
tolerance	float	Dimensionless bubble-point and log-fugacity tolerance.	1e-08
saturation_iterations	int	Maximum iterations for each bubble-point solve.	50
substitution_iterations	int	Batched sub-bubble flash iteration limits.	20
newton_iterations	int	Batched sub-bubble flash iteration limits.	20
raise_on_failure	bool	Raise :class:ConvergenceError if any boundary or flash fails. When false, failed-state viscosities are returned as NaN and diagnostics remain explicit.	True

Returns:

Type	Description
HeavyOilCSPProfile	Phase compositions, viscosities, and convergence diagnostics.

Raises:

Type	Description
ValueError:	If shapes, component order, or numerical controls are invalid.
InvalidStateError:	If state values or the feed are nonphysical.
ConvergenceError:	If a bubble point or sub-bubble flash fails and raise_on_failure=True.

Notes

The viscosity formulation is Lindeloff et al., "The Corresponding States Viscosity Model Applied to Heavy Oil Systems" (2004), Paper 2003-150, through Pedersen, Christensen, and Shaikh (2024), section 10.1.2, doi:10.1201/9780429457418. Phase equilibrium is a separate cubic-EoS calculation; its model and characterization parameter set must therefore be reported with any result.

Source code in src/torch_flash/transport/heavy_oil.py

def evaluate_heavy_oil_corresponding_states_profile(
    model: StateModel,
    temperature: Tensor,
    pressure: Tensor,
    feed_composition: Tensor,
    components: ComponentSet,
    *,
    third_csp: Tensor | float = 1.0,
    fourth_csp: Tensor | float = 1.0,
    initial_bubble_pressure: Tensor | float | None = None,
    tolerance: float = 1.0e-8,
    saturation_iterations: int = 50,
    substitution_iterations: int = 20,
    newton_iterations: int = 20,
    raise_on_failure: bool = True,
) -> HeavyOilCSPProfile:
    """Evaluate a phase-aware live-heavy-oil viscosity pressure profile.

    A liquid-vapor bubble point is solved once per unique temperature. States
    at or above that boundary retain the feed composition. All sub-bubble
    states at the same temperature are flashed simultaneously with
    :func:`torch_flash.batched_two_phase_flash`, and their equilibrium liquid
    compositions are passed in one batch to the heavy-oil CSP correlation.

    Parameters
    ----------
    model:
        Homogeneous-state fugacity model used for bubble points and flashes.
    temperature, pressure:
        One-dimensional state batches in K and Pa with matching shapes.
    feed_composition:
        One strictly positive feed mole-fraction vector.
    components:
        Transport critical properties and molar masses in the same order as
        the model and feed.
    third_csp, fourth_csp:
        Dimensionless Lindeloff calibration factors. Tensor inputs preserve
        gradients through the viscosity calculation.
    initial_bubble_pressure:
        Optional positive scalar pressure estimate in Pa reused at each
        temperature. Wilson initialization is used when omitted.
    tolerance:
        Dimensionless bubble-point and log-fugacity tolerance.
    saturation_iterations:
        Maximum iterations for each bubble-point solve.
    substitution_iterations, newton_iterations:
        Batched sub-bubble flash iteration limits.
    raise_on_failure:
        Raise :class:`ConvergenceError` if any boundary or flash fails. When
        false, failed-state viscosities are returned as NaN and diagnostics
        remain explicit.

    Returns
    -------
    HeavyOilCSPProfile
        Phase compositions, viscosities, and convergence diagnostics.

    Raises
    ------
    ValueError:
        If shapes, component order, or numerical controls are invalid.
    InvalidStateError:
        If state values or the feed are nonphysical.
    ConvergenceError:
        If a bubble point or sub-bubble flash fails and
        ``raise_on_failure=True``.

    Notes
    -----
    The viscosity formulation is Lindeloff et al., "The Corresponding States
    Viscosity Model Applied to Heavy Oil Systems" (2004), Paper 2003-150,
    through Pedersen, Christensen, and Shaikh (2024), section 10.1.2,
    doi:10.1201/9780429457418. Phase equilibrium is a separate cubic-EoS
    calculation; its model and characterization parameter set must therefore
    be reported with any result.
    """
    if (
        temperature.ndim != 1
        or pressure.ndim != 1
        or temperature.shape != pressure.shape
        or temperature.numel() == 0
    ):
        raise ValueError(
            "temperature and pressure must be nonempty matching one-dimensional batches"
        )
    if feed_composition.ndim != 1 or feed_composition.numel() != components.ncomponents:
        raise ValueError("feed composition must match the component set")
    if hasattr(model, "names") and tuple(model.names) != components.names:
        raise ValueError("phase model and transport component order must match")
    if (
        tolerance <= 0.0
        or saturation_iterations <= 0
        or substitution_iterations < 0
        or newton_iterations < 0
    ):
        raise ValueError("phase-equilibrium tolerances and iteration controls are invalid")
    invalid_state = (
        (~torch.isfinite(temperature))
        | (~torch.isfinite(pressure))
        | (temperature <= 0.0)
        | (pressure <= 0.0)
    )
    if bool(invalid_state.any()):
        raise InvalidStateError("heavy-oil profile temperatures and pressures must be positive")
    if not bool(
        torch.isfinite(feed_composition).all()
        & (feed_composition > 0.0).all()
        & torch.isfinite(feed_composition.sum())
    ):
        raise InvalidStateError("heavy-oil feed fractions must be finite and strictly positive")
    feed = normalize_composition(feed_composition)
    initial_pressure = (
        None
        if initial_bubble_pressure is None
        else torch.as_tensor(
            initial_bubble_pressure,
            dtype=temperature.dtype,
            device=temperature.device,
        )
    )
    if initial_pressure is not None and (
        initial_pressure.ndim
        or not bool(torch.isfinite(initial_pressure) & (initial_pressure > 0.0))
    ):
        raise ValueError("initial bubble pressure must be a positive finite scalar")

    state_count = temperature.numel()
    liquid_composition = feed.expand(state_count, -1).clone()
    bubble_pressure = torch.full_like(temperature, torch.nan)
    vapor_fraction = torch.zeros_like(temperature)
    residual_norm = torch.zeros_like(temperature)
    converged = torch.ones(state_count, dtype=torch.bool, device=temperature.device)

    for current_temperature in torch.unique(temperature):
        temperature_mask = temperature == current_temperature
        boundary = saturation_point(
            model,
            current_temperature,
            feed,
            "bubble",
            initial_pressure=initial_pressure,
            tolerance=tolerance,
            max_iterations=saturation_iterations,
        )
        bubble_pressure[temperature_mask] = boundary.pressure
        if not boundary.converged:
            converged[temperature_mask] = False
            residual_norm[temperature_mask] = boundary.residual_norm
            continue
        two_phase_mask = temperature_mask & (pressure < boundary.pressure)
        positions = torch.where(two_phase_mask)[0]
        if positions.numel() == 0:
            continue
        flash = batched_two_phase_flash(
            model,
            ChemicalState(
                temperature[positions],
                pressure[positions],
                feed,
            ),
            initial_k_values=boundary.k_values.expand(positions.numel(), -1),
            tolerance=tolerance,
            substitution_iterations=substitution_iterations,
            newton_iterations=newton_iterations,
        )
        liquid_composition[positions] = flash.liquid_composition
        vapor_fraction[positions] = flash.vapor_fraction
        residual_norm[positions] = flash.residual_norm
        converged[positions] = flash.converged

    if raise_on_failure and not bool(converged.all()):
        failed = torch.where(~converged)[0].detach().cpu().tolist()
        worst = float(residual_norm[~converged].detach().max())
        raise ConvergenceError(
            f"heavy-oil phase preparation failed at state indices {failed} "
            f"(maximum residual {worst:.3e})"
        )
    viscosity = heavy_oil_corresponding_states_viscosity(
        temperature,
        pressure,
        liquid_composition,
        components,
        phase="liquid",
        third_csp=third_csp,
        fourth_csp=fourth_csp,
    )
    viscosity = torch.where(converged, viscosity, torch.full_like(viscosity, torch.nan))
    return HeavyOilCSPProfile(
        viscosity,
        liquid_composition,
        bubble_pressure,
        vapor_fraction,
        converged,
        residual_norm,
    )

fit_heavy_oil_csp_factors

fit_heavy_oil_csp_factors(
    temperature: Tensor,
    pressure: Tensor,
    liquid_composition: Tensor,
    components: ComponentSet,
    observed_viscosity: Tensor,
    *,
    initial_factors: tuple[float, float] = (1.0, 1.0),
    factor_bounds: tuple[float, float] = (0.1, 5.0),
    optimizer: OptimizerKind = "lbfgs",
    learning_rate: float = 0.3,
    max_iterations: int = 200,
    tolerance: float = 1e-10,
    no_improvement_patience: int = 30,
) -> HeavyOilCSPCalibrationResult

Fit both heavy-oil CSP factors to all supplied states simultaneously.

The full-batch objective is the mean squared logarithmic viscosity ratio, so every experimental point contributes in every optimizer step without a unit-dependent scale. Physical factors are smoothly bounded with a sigmoid transform. PyTorch supplies optimizer gradients and the final local sensitivity matrix.

Parameters:

Name	Type	Description	Default
temperature	Tensor	State batches in K and Pa.	required
pressure	Tensor	State batches in K and Pa.	required
liquid_composition	Tensor	Feed or equilibrium liquid composition for every state, with components on the final axis.	required
components	ComponentSet	Critical properties and molar masses in matching order.	required
observed_viscosity	Tensor	Positive experimental dynamic viscosities in Pa s.	required
initial_factors	tuple[float, float]	Initial third and fourth CSP factors.	(1.0, 1.0)
factor_bounds	tuple[float, float]	Shared positive lower and upper bounds for both factors.	(0.1, 5.0)
optimizer	OptimizerKind	PyTorch optimizer identifier. Full-batch LBFGS with strong-Wolfe line search is the default for this two-parameter smooth problem.	'lbfgs'
learning_rate	float	Optimizer learning rate.	0.3
max_iterations	int	Maximum optimizer steps; the default is 200.	200
tolerance	float	Relative successive-loss stopping tolerance.	1e-10
no_improvement_patience	int	Early-stop patience; the default is 30 evaluated iterates.	30

Returns:

Type	Description
HeavyOilCSPCalibrationResult	Selected factors, fitted predictions, history, and identifiability diagnostics.

Raises:

Type	Description
ValueError:	If shapes, bounds, factors, or optimizer controls are invalid.
InvalidStateError:	If observations are nonpositive or nonfinite.

Notes

The fitted factors are the multipliers in Pedersen, Christensen, and Shaikh (2024), Eqs. 10.34-10.35, doi:10.1201/9780429457418. Lindeloff et al. (2004), Paper 2003-150, likewise fit both factors to all three Oil 5 temperature series simultaneously. This API does not refit the phase-equilibrium model.

Source code in src/torch_flash/transport/heavy_oil.py

def fit_heavy_oil_csp_factors(
    temperature: Tensor,
    pressure: Tensor,
    liquid_composition: Tensor,
    components: ComponentSet,
    observed_viscosity: Tensor,
    *,
    initial_factors: tuple[float, float] = (1.0, 1.0),
    factor_bounds: tuple[float, float] = (0.1, 5.0),
    optimizer: OptimizerKind = "lbfgs",
    learning_rate: float = 0.3,
    max_iterations: int = 200,
    tolerance: float = 1.0e-10,
    no_improvement_patience: int = 30,
) -> HeavyOilCSPCalibrationResult:
    """Fit both heavy-oil CSP factors to all supplied states simultaneously.

    The full-batch objective is the mean squared logarithmic viscosity ratio,
    so every experimental point contributes in every optimizer step without
    a unit-dependent scale. Physical factors are smoothly bounded with a
    sigmoid transform. PyTorch supplies optimizer gradients and the final
    local sensitivity matrix.

    Parameters
    ----------
    temperature, pressure:
        State batches in K and Pa.
    liquid_composition:
        Feed or equilibrium liquid composition for every state, with
        components on the final axis.
    components:
        Critical properties and molar masses in matching order.
    observed_viscosity:
        Positive experimental dynamic viscosities in Pa s.
    initial_factors:
        Initial third and fourth CSP factors.
    factor_bounds:
        Shared positive lower and upper bounds for both factors.
    optimizer:
        PyTorch optimizer identifier. Full-batch LBFGS with strong-Wolfe line
        search is the default for this two-parameter smooth problem.
    learning_rate:
        Optimizer learning rate.
    max_iterations:
        Maximum optimizer steps; the default is 200.
    tolerance:
        Relative successive-loss stopping tolerance.
    no_improvement_patience:
        Early-stop patience; the default is 30 evaluated iterates.

    Returns
    -------
    HeavyOilCSPCalibrationResult
        Selected factors, fitted predictions, history, and identifiability
        diagnostics.

    Raises
    ------
    ValueError:
        If shapes, bounds, factors, or optimizer controls are invalid.
    InvalidStateError:
        If observations are nonpositive or nonfinite.

    Notes
    -----
    The fitted factors are the multipliers in Pedersen, Christensen, and
    Shaikh (2024), Eqs. 10.34-10.35,
    doi:10.1201/9780429457418. Lindeloff et al. (2004), Paper 2003-150,
    likewise fit both factors to all three Oil 5 temperature series
    simultaneously. This API does not refit the phase-equilibrium model.
    """
    if (
        temperature.ndim != 1
        or pressure.shape != temperature.shape
        or observed_viscosity.shape != temperature.shape
        or liquid_composition.shape != (temperature.numel(), components.ncomponents)
    ):
        raise ValueError("heavy-oil calibration state, observation, and composition shapes differ")
    if not bool(torch.isfinite(observed_viscosity).all() & (observed_viscosity > 0.0).all()):
        raise InvalidStateError("observed viscosities must be finite and positive")
    lower, upper = factor_bounds
    if (
        not 0.0 < lower < upper
        or len(initial_factors) != 2
        or any(not lower < value < upper for value in initial_factors)
    ):
        raise ValueError("initial CSP factors must lie strictly inside positive ordered bounds")
    lower_tensor = temperature.new_tensor(lower)
    span = temperature.new_tensor(upper - lower)
    initial = temperature.new_tensor(initial_factors)
    initial_unit = (initial - lower_tensor) / span
    unconstrained = nn.Parameter(torch.logit(initial_unit))

    def physical_factors(raw: Tensor) -> Tensor:
        return lower_tensor + span * torch.sigmoid(raw)

    def closure() -> Tensor:
        factors = physical_factors(unconstrained)
        prediction = heavy_oil_corresponding_states_viscosity(
            temperature,
            pressure,
            liquid_composition,
            components,
            phase="liquid",
            third_csp=factors[0],
            fourth_csp=factors[1],
        )
        return torch.mean(torch.log(prediction / observed_viscosity).square())

    fit = fit_parameters(
        (unconstrained,),
        closure,
        optimizer=optimizer,
        learning_rate=learning_rate,
        max_iterations=max_iterations,
        tolerance=tolerance,
        no_improvement_patience=no_improvement_patience,
    )
    factors = physical_factors(unconstrained).detach()
    prediction = heavy_oil_corresponding_states_viscosity(
        temperature,
        pressure,
        liquid_composition,
        components,
        phase="liquid",
        third_csp=factors[0],
        fourth_csp=factors[1],
    ).detach()

    def log_prediction(current_factors: Tensor) -> Tensor:
        return torch.log(
            heavy_oil_corresponding_states_viscosity(
                temperature,
                pressure,
                liquid_composition,
                components,
                phase="liquid",
                third_csp=current_factors[0],
                fourth_csp=current_factors[1],
            )
        )

    sensitivity = torch.func.jacrev(log_prediction)(factors)
    singular_values = torch.linalg.svdvals(sensitivity)
    rank_tolerance = (
        max(sensitivity.shape) * torch.finfo(sensitivity.dtype).eps * singular_values.max()
    )
    rank = int(torch.sum(singular_values > rank_tolerance))
    condition_number = (
        float(torch.inf) if rank < 2 else float((singular_values[0] / singular_values[-1]).detach())
    )
    return HeavyOilCSPCalibrationResult(
        factors[0],
        factors[1],
        prediction,
        fit,
        sensitivity,
        singular_values,
        rank,
        condition_number,
    )

brock_bird_surface_tension

brock_bird_surface_tension(
    temperature: Tensor,
    critical_temperature: Tensor,
    critical_pressure: Tensor,
    normal_boiling_temperature: Tensor,
) -> Tensor

Evaluate Brock-Bird pure nonpolar-fluid surface tension.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperatures in K.	required
critical_temperature	Tensor	Temperatures in K.	required
normal_boiling_temperature	Tensor	Temperatures in K.	required
critical_pressure	Tensor	Critical pressure in Pa.	required

Returns:

Type	Description
Tensor	Surface tension in N/m.

Raises:

Type	Description
InvalidStateError	If the state is not below the critical temperature or the defining properties are nonphysical.

Notes

Implements Brock and Bird, "Surface Tension and the Principle of Corresponding States," AIChE Journal 1 (1955), 174-177, through Pedersen et al. (2024), Eqs. 10.77-10.79. The published numerical equation yields dyn/cm from K and atm; the result is converted to N/m.

Source code in src/torch_flash/transport/interfacial_tension.py

def brock_bird_surface_tension(
    temperature: Tensor,
    critical_temperature: Tensor,
    critical_pressure: Tensor,
    normal_boiling_temperature: Tensor,
) -> Tensor:
    """Evaluate Brock-Bird pure nonpolar-fluid surface tension.

    Parameters
    ----------
    temperature, critical_temperature, normal_boiling_temperature
        Temperatures in K.
    critical_pressure
        Critical pressure in Pa.

    Returns
    -------
    Tensor
        Surface tension in N/m.

    Raises
    ------
    InvalidStateError
        If the state is not below the critical temperature or the defining
        properties are nonphysical.

    Notes
    -----
    Implements Brock and Bird, "Surface Tension and the Principle of
    Corresponding States," *AIChE Journal* 1 (1955), 174-177, through
    Pedersen et al. (2024), Eqs. 10.77-10.79. The published numerical
    equation yields dyn/cm from K and atm; the result is converted to N/m.
    """
    temperature, critical, pressure, boiling = torch.broadcast_tensors(
        temperature,
        critical_temperature,
        critical_pressure,
        normal_boiling_temperature,
    )
    if bool(
        ((~torch.isfinite(temperature)) | (temperature <= 0.0) | (temperature >= critical)).any()
    ):
        raise InvalidStateError("Brock-Bird requires a positive subcritical temperature")
    alpha = riedel_parameter(boiling, critical, pressure)
    amplitude = (pressure / 101_325.0).pow(2.0 / 3.0) * critical.pow(1.0 / 3.0)
    amplitude = amplitude * (0.133 * alpha - 0.281)
    tension_dyn_cm = amplitude * (1.0 - temperature / critical).pow(11.0 / 9.0)
    if bool((~torch.isfinite(tension_dyn_cm) | (tension_dyn_cm < 0.0)).any()):
        raise InvalidStateError("Brock-Bird produced a negative surface tension")
    result: Tensor = 1.0e-3 * tension_dyn_cm
    return result

lee_chien_interfacial_tension

lee_chien_interfacial_tension(
    liquid_molar_density: Tensor,
    vapor_molar_density: Tensor,
    liquid_composition: Tensor,
    vapor_composition: Tensor,
    components: ComponentSet,
    riedel_parameters: Tensor,
    b_parameters: Tensor,
) -> Tensor

Evaluate the Lee-Chien phase-parachor interfacial tension.

Parameters:

Name	Type	Description	Default
liquid_molar_density	Tensor	Coexisting molar densities in mol/m3.	required
vapor_molar_density	Tensor	Coexisting molar densities in mol/m3.	required
liquid_composition	Tensor	Phase mole fractions with components on the final axis.	required
vapor_composition	Tensor	Phase mole fractions with components on the final axis.	required
components	ComponentSet	Critical temperatures, pressures, and molar volumes.	required
riedel_parameters	Tensor	Dimensionless per-component Riedel parameters.	required
b_parameters	Tensor	Per-component Lee-Chien B parameters.	required

Returns:

Type	Description
Tensor	Interfacial tension in N/m.

Notes

Implements Lee and Chien, SPE-12643-MS (1984), through Pedersen et al. (2024), Eqs. 10.86-10.97. The bundled B table is deliberately limited to the eight components reproduced in Pedersen Table 10.20; additional components require explicit values. C7+ whole-fraction construction is a characterization step and is not inferred silently.

Source code in src/torch_flash/transport/interfacial_tension.py

def lee_chien_interfacial_tension(
    liquid_molar_density: Tensor,
    vapor_molar_density: Tensor,
    liquid_composition: Tensor,
    vapor_composition: Tensor,
    components: ComponentSet,
    riedel_parameters: Tensor,
    b_parameters: Tensor,
) -> Tensor:
    """Evaluate the Lee-Chien phase-parachor interfacial tension.

    Parameters
    ----------
    liquid_molar_density, vapor_molar_density
        Coexisting molar densities in mol/m3.
    liquid_composition, vapor_composition
        Phase mole fractions with components on the final axis.
    components
        Critical temperatures, pressures, and molar volumes.
    riedel_parameters
        Dimensionless per-component Riedel parameters.
    b_parameters
        Per-component Lee-Chien ``B`` parameters.

    Returns
    -------
    Tensor
        Interfacial tension in N/m.

    Notes
    -----
    Implements Lee and Chien, SPE-12643-MS (1984), through Pedersen et al.
    (2024), Eqs. 10.86-10.97. The bundled ``B`` table is deliberately limited
    to the eight components reproduced in Pedersen Table 10.20; additional
    components require explicit values. C7+ whole-fraction construction is a
    characterization step and is not inferred silently.
    """
    if liquid_composition.shape != vapor_composition.shape:
        raise ValueError("liquid and vapor compositions must have the same shape")
    if liquid_composition.ndim < 1 or liquid_composition.shape[-1] != components.ncomponents:
        raise ValueError("composition and component set sizes must match")
    expected = components.critical_temperature.shape
    if riedel_parameters.shape != expected or b_parameters.shape != expected:
        raise ValueError("Riedel and B parameters must have one value per component")
    if components.critical_volume is None:
        raise ValueError("Lee-Chien requires component critical molar volumes")
    critical_volume = components.critical_volume
    if bool(
        (
            (~torch.isfinite(riedel_parameters))
            | (~torch.isfinite(b_parameters))
            | (b_parameters <= 0.0)
        ).any()
    ):
        raise ValueError("Lee-Chien parameters must be finite and B positive")
    liquid = normalize_composition(liquid_composition)
    vapor = normalize_composition(vapor_composition)

    def phase_parachor(composition: Tensor) -> Tensor:
        pressure_atm = torch.sum(
            composition * (components.critical_pressure / 101_325.0),
            dim=-1,
        )
        critical_temperature = torch.sum(
            composition * components.critical_temperature,
            dim=-1,
        )
        alpha = torch.sum(composition * riedel_parameters, dim=-1)
        critical_volume_cm3 = torch.sum(
            composition * critical_volume * 1.0e6,
            dim=-1,
        )
        b = torch.sum(composition * b_parameters, dim=-1)
        amplitude = (
            pressure_atm.pow(2.0 / 3.0)
            * critical_temperature.pow(1.0 / 3.0)
            * (0.133 * alpha - 0.281)
        )
        return amplitude.pow(0.25) * critical_volume_cm3 / b

    liquid_parachor = phase_parachor(liquid)
    vapor_parachor = phase_parachor(vapor)
    rho_l, rho_v = torch.broadcast_tensors(liquid_molar_density, vapor_molar_density)
    if bool(
        (
            (~torch.isfinite(rho_l)) | (~torch.isfinite(rho_v)) | (rho_l <= rho_v) | (rho_v < 0.0)
        ).any()
    ):
        raise InvalidStateError("Lee-Chien requires rho_liquid > rho_vapor >= 0")
    contrast = (rho_l / 1.0e6 * liquid_parachor - rho_v / 1.0e6 * vapor_parachor).abs()
    result: Tensor = 1.0e-3 * contrast.pow(4)
    if bool((~torch.isfinite(result) | (result < 0.0)).any()):
        raise InvalidStateError("Lee-Chien produced a non-finite interfacial tension")
    return result

parachor_from_molar_mass

parachor_from_molar_mass(molar_mass: Tensor) -> Tensor

Estimate a C7+ petroleum-fraction parachor from molar mass.

Parameters:

Name	Type	Description	Default
molar_mass	Tensor	Fraction molar mass in kg/mol.	required

Returns:

Type	Description
Tensor	Parachor in the Weinaug-Katz (dyn/cm)**0.25 cm3/mol convention.

Notes

Implements Pedersen et al. (2024), Eq. 10.85. This is a petroleum pseudo-component correlation, not a replacement for published pure-component parachors.

Source code in src/torch_flash/transport/interfacial_tension.py

def parachor_from_molar_mass(molar_mass: Tensor) -> Tensor:
    """Estimate a C7+ petroleum-fraction parachor from molar mass.

    Parameters
    ----------
    molar_mass
        Fraction molar mass in kg/mol.

    Returns
    -------
    Tensor
        Parachor in the Weinaug-Katz ``(dyn/cm)**0.25 cm3/mol`` convention.

    Notes
    -----
    Implements Pedersen et al. (2024), Eq. 10.85. This is a petroleum
    pseudo-component correlation, not a replacement for published
    pure-component parachors.
    """
    if bool(((~torch.isfinite(molar_mass)) | (molar_mass <= 0.0)).any()):
        raise InvalidStateError("molar mass must be finite and positive")
    return 59.3 + 2.34 * (1000.0 * molar_mass)

published_lee_chien_b

published_lee_chien_b(
    names: Sequence[str], *, like: Tensor
) -> Tensor

Return Lee-Chien B values reproduced in Pedersen Table 10.20.

Parameters:

Name	Type	Description	Default
names	Sequence[str]	Canonical component names in the desired output order.	required
like	Tensor	Tensor supplying output dtype and device.	required

Returns:

Type	Description
Tensor	One dimensionless Lee-Chien B coefficient per name.

Raises:

Type	Description
KeyError	If a component lies outside the eight-component published table. Explicit user-supplied values are required for additional compounds.

Notes

The coefficients are reproduced from Pedersen, Christensen, and Shaikh, Phase Behavior of Petroleum Reservoir Fluids, 3rd ed. (2024), Table 10.20, doi:10.1201/9780429457418. The defining correlation is Lee and Chien, SPE-12643-MS (1984), doi:10.2118/12643-MS.

Source code in src/torch_flash/transport/interfacial_tension.py

def published_lee_chien_b(
    names: Sequence[str],
    *,
    like: Tensor,
) -> Tensor:
    """Return Lee-Chien ``B`` values reproduced in Pedersen Table 10.20.

    Parameters
    ----------
    names
        Canonical component names in the desired output order.
    like
        Tensor supplying output dtype and device.

    Returns
    -------
    Tensor
        One dimensionless Lee-Chien ``B`` coefficient per name.

    Raises
    ------
    KeyError
        If a component lies outside the eight-component published table.
        Explicit user-supplied values are required for additional compounds.

    Notes
    -----
    The coefficients are reproduced from Pedersen, Christensen, and Shaikh,
    *Phase Behavior of Petroleum Reservoir Fluids*, 3rd ed. (2024),
    Table 10.20, doi:10.1201/9780429457418. The defining correlation is Lee
    and Chien, SPE-12643-MS (1984), doi:10.2118/12643-MS.
    """
    values = cast(Mapping[str, float], _INTERFACIAL["lee_chien_b"])
    missing = [name for name in names if name not in values]
    if missing:
        raise KeyError(f"no published Lee-Chien B value for: {', '.join(missing)}")
    return like.new_tensor([values[name] for name in names])

published_parachors

published_parachors(
    names: Sequence[str], *, like: Tensor
) -> Tensor

Return the Pedersen Table 10.19 parachors for supported components.

Parameters:

Name	Type	Description	Default
names	Sequence[str]	Canonical component names in the desired output order.	required
like	Tensor	Tensor supplying output dtype and device.	required

Returns:

Type	Description
Tensor	One published component parachor per name, in conventional (dyn/cm)**0.25 cm3/mol units.

Raises:

Type	Description
KeyError	If any component is outside the published table. Callers may instead provide explicit parachors or use :func:parachor_from_molar_mass for characterized C7+ fractions.

Notes

Values are reproduced from Pedersen, Christensen, and Shaikh, Phase Behavior of Petroleum Reservoir Fluids, 3rd ed. (2024), Table 10.19, doi:10.1201/9780429457418.

Source code in src/torch_flash/transport/interfacial_tension.py

def published_parachors(
    names: Sequence[str],
    *,
    like: Tensor,
) -> Tensor:
    """Return the Pedersen Table 10.19 parachors for supported components.

    Parameters
    ----------
    names
        Canonical component names in the desired output order.
    like
        Tensor supplying output dtype and device.

    Returns
    -------
    Tensor
        One published component parachor per name, in conventional
        ``(dyn/cm)**0.25 cm3/mol`` units.

    Raises
    ------
    KeyError
        If any component is outside the published table. Callers may instead
        provide explicit parachors or use :func:`parachor_from_molar_mass` for
        characterized C7+ fractions.

    Notes
    -----
    Values are reproduced from Pedersen, Christensen, and Shaikh, *Phase
    Behavior of Petroleum Reservoir Fluids*, 3rd ed. (2024), Table 10.19,
    doi:10.1201/9780429457418.
    """
    values = cast(Mapping[str, float], _INTERFACIAL["parachor"])
    missing = [name for name in names if name not in values]
    if missing:
        raise KeyError(f"no published parachor for: {', '.join(missing)}")
    return like.new_tensor([values[name] for name in names])

riedel_parameter

riedel_parameter(
    normal_boiling_temperature: Tensor,
    critical_temperature: Tensor,
    critical_pressure: Tensor,
) -> Tensor

Return the dimensionless Riedel parameter.

Parameters:

Name	Type	Description	Default
normal_boiling_temperature	Tensor	Normal boiling and critical temperatures in K.	required
critical_temperature	Tensor	Normal boiling and critical temperatures in K.	required
critical_pressure	Tensor	Critical pressure in Pa.	required

Returns:

Type	Description
Tensor	Dimensionless Riedel parameter.

Notes

Implements Pedersen et al. (2024), Eq. 10.79, from Riedel (1954). The numerical form uses critical pressure in atm internally.

Source code in src/torch_flash/transport/interfacial_tension.py

def riedel_parameter(
    normal_boiling_temperature: Tensor,
    critical_temperature: Tensor,
    critical_pressure: Tensor,
) -> Tensor:
    """Return the dimensionless Riedel parameter.

    Parameters
    ----------
    normal_boiling_temperature, critical_temperature
        Normal boiling and critical temperatures in K.
    critical_pressure
        Critical pressure in Pa.

    Returns
    -------
    Tensor
        Dimensionless Riedel parameter.

    Notes
    -----
    Implements Pedersen et al. (2024), Eq. 10.79, from Riedel (1954). The
    numerical form uses critical pressure in atm internally.
    """
    boiling, critical, pressure = torch.broadcast_tensors(
        normal_boiling_temperature,
        critical_temperature,
        critical_pressure,
    )
    if bool(
        (
            (~torch.isfinite(boiling))
            | (~torch.isfinite(critical))
            | (~torch.isfinite(pressure))
            | (boiling <= 0.0)
            | (critical <= 0.0)
            | (boiling >= critical)
            | (pressure <= 0.0)
        ).any()
    ):
        raise InvalidStateError("Riedel inputs require 0 < Tb < Tc and positive Pc")
    reduced_boiling = boiling / critical
    result: Tensor = 0.9076 * (
        1.0 + reduced_boiling * torch.log(pressure / 101_325.0) / (1.0 - reduced_boiling)
    )
    return result

weinaug_katz_interfacial_tension

weinaug_katz_interfacial_tension(
    liquid_molar_density: Tensor,
    vapor_molar_density: Tensor,
    liquid_composition: Tensor,
    vapor_composition: Tensor,
    parachor: Tensor,
    *,
    liquid_mass_density: Tensor | None = None,
    vapor_mass_density: Tensor | None = None,
    danesh_exponent: bool = False,
) -> Tensor

Evaluate Weinaug-Katz gas-oil interfacial tension.

Parameters:

Name	Type	Description	Default
liquid_molar_density	Tensor	Coexisting phase molar densities in mol/m3.	required
vapor_molar_density	Tensor	Coexisting phase molar densities in mol/m3.	required
liquid_composition	Tensor	Phase mole fractions with components on the final axis.	required
vapor_composition	Tensor	Phase mole fractions with components on the final axis.	required
parachor	Tensor	Component parachors in the conventional (dyn/cm)**0.25 cm3/mol units.	required
liquid_mass_density	Tensor | None	Optional phase mass densities in kg/m3 required by the Danesh gas-condensate exponent.	None
vapor_mass_density	Tensor | None	Optional phase mass densities in kg/m3 required by the Danesh gas-condensate exponent.	None
danesh_exponent	bool	Use Pedersen Eq. 10.98 instead of the fixed one-quarter exponent.	False

Returns:

Type	Description
Tensor	Interfacial tension in N/m.

Raises:

Type	Description
ValueError	If component dimensions disagree or Danesh densities are absent.
InvalidStateError	If phase densities are nonphysical.

Notes

Implements the Macleod-Sugden/Weinaug-Katz relation as reproduced in Pedersen et al. (2024), Eqs. 10.83-10.85. The optional Danesh et al. modification follows Eq. 10.98. This operation requires an already solved pair of coexisting phases and does not identify or flash phases.

Source code in src/torch_flash/transport/interfacial_tension.py

def weinaug_katz_interfacial_tension(
    liquid_molar_density: Tensor,
    vapor_molar_density: Tensor,
    liquid_composition: Tensor,
    vapor_composition: Tensor,
    parachor: Tensor,
    *,
    liquid_mass_density: Tensor | None = None,
    vapor_mass_density: Tensor | None = None,
    danesh_exponent: bool = False,
) -> Tensor:
    """Evaluate Weinaug-Katz gas-oil interfacial tension.

    Parameters
    ----------
    liquid_molar_density, vapor_molar_density
        Coexisting phase molar densities in mol/m3.
    liquid_composition, vapor_composition
        Phase mole fractions with components on the final axis.
    parachor
        Component parachors in the conventional
        ``(dyn/cm)**0.25 cm3/mol`` units.
    liquid_mass_density, vapor_mass_density
        Optional phase mass densities in kg/m3 required by the Danesh
        gas-condensate exponent.
    danesh_exponent
        Use Pedersen Eq. 10.98 instead of the fixed one-quarter exponent.

    Returns
    -------
    Tensor
        Interfacial tension in N/m.

    Raises
    ------
    ValueError
        If component dimensions disagree or Danesh densities are absent.
    InvalidStateError
        If phase densities are nonphysical.

    Notes
    -----
    Implements the Macleod-Sugden/Weinaug-Katz relation as reproduced in
    Pedersen et al. (2024), Eqs. 10.83-10.85. The optional Danesh et al.
    modification follows Eq. 10.98. This operation requires an already solved
    pair of coexisting phases and does not identify or flash phases.
    """
    if liquid_composition.shape != vapor_composition.shape:
        raise ValueError("liquid and vapor compositions must have the same shape")
    if liquid_composition.ndim < 1 or parachor.shape != liquid_composition.shape[-1:]:
        raise ValueError("parachor must have one value per composition component")
    liquid = normalize_composition(liquid_composition)
    vapor = normalize_composition(vapor_composition)
    rho_l, rho_v = torch.broadcast_tensors(liquid_molar_density, vapor_molar_density)
    if bool(
        (
            (~torch.isfinite(rho_l))
            | (~torch.isfinite(rho_v))
            | (rho_l <= 0.0)
            | (rho_v < 0.0)
            | (rho_l <= rho_v)
        ).any()
    ):
        raise InvalidStateError("interfacial tension requires rho_liquid > rho_vapor >= 0")
    density_liquid_mol_cm3 = rho_l / 1.0e6
    density_vapor_mol_cm3 = rho_v / 1.0e6
    contrast = torch.sum(
        parachor
        * (density_liquid_mol_cm3[..., None] * liquid - density_vapor_mol_cm3[..., None] * vapor),
        dim=-1,
    ).abs()
    exponent = contrast.new_tensor(0.25)
    if danesh_exponent:
        if liquid_mass_density is None or vapor_mass_density is None:
            raise ValueError("Danesh exponent requires liquid and vapor mass densities")
        mass_liquid, mass_vapor = torch.broadcast_tensors(
            liquid_mass_density,
            vapor_mass_density,
        )
        if bool(
            (
                (~torch.isfinite(mass_liquid))
                | (~torch.isfinite(mass_vapor))
                | (mass_liquid <= mass_vapor)
                | (mass_vapor < 0.0)
            ).any()
        ):
            raise InvalidStateError("Danesh exponent requires rho_liquid > rho_vapor >= 0")
        density_difference_g_cm3 = (mass_liquid - mass_vapor) / 1000.0
        exponent = 1.0 / (3.583 + 0.16 * density_difference_g_cm3)
    tension_dyn_cm = contrast.pow(1.0 / exponent)
    return 1.0e-3 * tension_dyn_cm

corresponding_states_thermal_conductivity

corresponding_states_thermal_conductivity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    components: ComponentSet,
    ideal_gas: HeatCapacityModel,
    methane_ideal_gas: HeatCapacityModel,
    *,
    phase: Literal["liquid", "vapor"] = "vapor",
) -> Tensor

Evaluate the Christensen-Fredenslund mixture thermal conductivity.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Positive state variables in K and Pa.	required
pressure	Tensor	Positive state variables in K and Pa.	required
composition	Tensor	Mole fractions with components on the final axis.	required
components	ComponentSet	Critical properties and molar masses matching composition.	required
ideal_gas	HeatCapacityModel	Component ideal-gas heat-capacity model for the mixture.	required
methane_ideal_gas	HeatCapacityModel	One-component methane ideal-gas heat-capacity model.	required
phase	Literal['liquid', 'vapor']	Methane BWR reference-density branch.	'vapor'

Returns:

Type	Description
Tensor	Mixture thermal conductivity in W/(m K).

Raises:

Type	Description
ValueError	If composition or heat-capacity component dimensions disagree.
InvalidStateError	If the model produces a nonpositive result.

Notes

Implements Christensen and Fredenslund, "A corresponding states model for the thermal conductivity of gases and liquids," Chemical Engineering Science 35 (1980), 871-875, and Pedersen et al. (2024), Eqs. 10.59-10.74. As specified for mixtures, the pure-fluid critical enhancement is omitted. The model requires homogeneous phase identity; it does not perform a flash.

Source code in src/torch_flash/transport/thermal_conductivity.py

def corresponding_states_thermal_conductivity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    components: ComponentSet,
    ideal_gas: HeatCapacityModel,
    methane_ideal_gas: HeatCapacityModel,
    *,
    phase: Literal["liquid", "vapor"] = "vapor",
) -> Tensor:
    """Evaluate the Christensen-Fredenslund mixture thermal conductivity.

    Parameters
    ----------
    temperature, pressure
        Positive state variables in K and Pa.
    composition
        Mole fractions with components on the final axis.
    components
        Critical properties and molar masses matching ``composition``.
    ideal_gas
        Component ideal-gas heat-capacity model for the mixture.
    methane_ideal_gas
        One-component methane ideal-gas heat-capacity model.
    phase
        Methane BWR reference-density branch.

    Returns
    -------
    Tensor
        Mixture thermal conductivity in W/(m K).

    Raises
    ------
    ValueError
        If composition or heat-capacity component dimensions disagree.
    InvalidStateError
        If the model produces a nonpositive result.

    Notes
    -----
    Implements Christensen and Fredenslund, "A corresponding states model for
    the thermal conductivity of gases and liquids," *Chemical Engineering
    Science* 35 (1980), 871-875, and Pedersen et al. (2024),
    Eqs. 10.59-10.74. As specified for mixtures, the pure-fluid critical
    enhancement is omitted. The model requires homogeneous phase identity;
    it does not perform a flash.
    """
    if composition.ndim < 1 or composition.shape[-1] != components.ncomponents:
        raise ValueError("composition and component set sizes must match")
    z = normalize_composition(composition)
    mixture_tc, mixture_pc = _mixture_pseudocritical(components, z)
    mapped_temperature = temperature * _METHANE_TC / mixture_tc
    mapped_pressure = pressure * _METHANE_PC / mixture_pc
    initial_density = methane_bwr_density(mapped_temperature, mapped_pressure, phase=phase)
    reduced_density = (
        initial_density * _METHANE_MOLAR_MASS_G / 1000.0 / _METHANE_CRITICAL_MASS_DENSITY
    )
    molecular_weight_g = 1000.0 * components.molar_mass
    alpha_components = 1.0 + 0.0006004 * reduced_density[..., None].pow(
        2.043
    ) * molecular_weight_g.pow(1.086)
    weights = z[..., :, None] * z[..., None, :]
    alpha_mixture = torch.sum(
        weights * torch.sqrt(alpha_components[..., :, None] * alpha_components[..., None, :]),
        dim=(-2, -1),
    )
    alpha_methane = 1.0 + 0.0006004 * reduced_density.pow(2.043) * temperature.new_tensor(
        _METHANE_MOLAR_MASS_G
    ).pow(1.086)
    reference_temperature = temperature * _METHANE_TC * alpha_methane / (mixture_tc * alpha_mixture)
    reference_pressure = pressure * _METHANE_PC * alpha_methane / (mixture_pc * alpha_mixture)
    reference_density = methane_bwr_density(
        reference_temperature,
        reference_pressure,
        phase=phase,
    )
    reference_conductivity = methane_thermal_conductivity(
        reference_temperature,
        reference_density,
        include_critical_enhancement=False,
    )
    effective_mass = _effective_mixture_molar_mass(z, components, mixture_tc, mixture_pc)
    mixture_component_heat_capacity = ideal_gas.heat_capacity(temperature)
    if (
        mixture_component_heat_capacity.ndim < 1
        or mixture_component_heat_capacity.shape[-1] != components.ncomponents
    ):
        raise ValueError("mixture heat-capacity model must return one value per component")
    methane_component_heat_capacity = methane_ideal_gas.heat_capacity(reference_temperature)
    if methane_component_heat_capacity.ndim < 1 or methane_component_heat_capacity.shape[-1] != 1:
        raise ValueError("methane heat-capacity model must return one component value")
    mixture_heat_capacity = torch.sum(z * mixture_component_heat_capacity, dim=-1)
    methane_heat_capacity = methane_component_heat_capacity.squeeze(-1)
    mixture_dilute_viscosity = corresponding_states_viscosity(
        temperature,
        temperature.new_tensor(101_325.0),
        z,
        components,
        phase="vapor",
    )
    methane_dilute_viscosity = methane_viscosity(
        reference_temperature,
        methane_bwr_density(
            reference_temperature,
            reference_temperature.new_tensor(101_325.0),
            phase="vapor",
        ),
    )
    mixture_internal = _internal_energy_conductivity(
        mixture_dilute_viscosity,
        mixture_heat_capacity,
        effective_mass,
        reduced_density,
    )
    reference_internal = _internal_energy_conductivity(
        methane_dilute_viscosity,
        methane_heat_capacity,
        reference_temperature.new_tensor(_METHANE_MOLAR_MASS_G / 2.0),
        reduced_density,
    )
    scale = (
        (mixture_pc / _METHANE_PC).pow(2.0 / 3.0)
        / (
            (mixture_tc / _METHANE_TC).pow(1.0 / 6.0)
            * (effective_mass / (_METHANE_MOLAR_MASS_G / 2.0)).sqrt()
        )
        * alpha_mixture
        / alpha_methane
    )
    conductivity = scale * (reference_conductivity - reference_internal) + mixture_internal
    if bool((~torch.isfinite(conductivity) | (conductivity <= 0.0)).any()):
        raise InvalidStateError("corresponding-states model produced non-positive conductivity")
    return conductivity

methane_critical_thermal_conductivity_enhancement

methane_critical_thermal_conductivity_enhancement(
    temperature: Tensor, density_mol_l: Tensor
) -> Tensor

Return the methane critical thermal-conductivity enhancement.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Methane temperature in K.	required
density_mol_l	Tensor	Methane molar density in mol/L.	required

Returns:

Type	Description
Tensor	Critical enhancement in W/(m K).

Raises:

Type	Description
InvalidStateError	If the supplied state is nonpositive or has nonpositive isothermal compressibility.

Notes

Implements Hanley, McCarty, and Haynes, Cryogenics 15 (1975), Appendix C, Eqs. 10-17, doi:10.1016/0011-2275(75)90010-7. The Vicentini-Missoni scaling form replaces the BWR compressibility only inside |T-Tc|/Tc <= 0.025 and |rho-rhoc|/rhoc <= 0.25. The correlation is singular at the mathematical critical point.

Source code in src/torch_flash/transport/thermal_conductivity.py

def methane_critical_thermal_conductivity_enhancement(
    temperature: Tensor,
    density_mol_l: Tensor,
) -> Tensor:
    """Return the methane critical thermal-conductivity enhancement.

    Parameters
    ----------
    temperature
        Methane temperature in K.
    density_mol_l
        Methane molar density in mol/L.

    Returns
    -------
    Tensor
        Critical enhancement in W/(m K).

    Raises
    ------
    InvalidStateError
        If the supplied state is nonpositive or has nonpositive isothermal
        compressibility.

    Notes
    -----
    Implements Hanley, McCarty, and Haynes, *Cryogenics* 15 (1975),
    Appendix C, Eqs. 10-17, doi:10.1016/0011-2275(75)90010-7. The
    Vicentini-Missoni scaling form replaces the BWR compressibility only
    inside ``|T-Tc|/Tc <= 0.025`` and ``|rho-rhoc|/rhoc <= 0.25``. The
    correlation is singular at the mathematical critical point.
    """
    temperature, density = torch.broadcast_tensors(temperature, density_mol_l)
    if bool(
        (
            (~torch.isfinite(temperature))
            | (~torch.isfinite(density))
            | (temperature <= 0.0)
            | (density <= 0.0)
        ).any()
    ):
        raise InvalidStateError("methane critical enhancement requires positive finite state data")
    density_mol_m3 = 1000.0 * density
    derivative = _methane_bwr_density_derivative(temperature, density)
    pressure_density_derivative = 101_325.0 * derivative / 1000.0
    direct_compressibility = 1.0 / (density_mol_m3 * pressure_density_derivative)
    compressibility = _scaled_isothermal_compressibility(
        temperature,
        density,
        direct_compressibility,
    )
    if bool((~torch.isfinite(compressibility) | (compressibility <= 0.0)).any()):
        raise InvalidStateError("methane BWR state has non-positive isothermal compressibility")
    mass_density_g_cm3 = density * _METHANE_MOLAR_MASS_G / 1000.0
    mass_density_kg_m3 = 1000.0 * mass_density_g_cm3
    viscosity = methane_viscosity(temperature, density)
    pressure_temperature = _pressure_temperature_derivative(temperature, density)
    correlation_length = 1.465 * torch.sqrt(mass_density_g_cm3 / temperature) * 1.0e-8
    prefactor = torch.sqrt(
        (_METHANE_MOLAR_MASS_G / 1000.0)
        / (mass_density_kg_m3 * _AVOGADRO_CONSTANT * _BOLTZMANN_CONSTANT * temperature)
    )
    reduced_temperature_distance = torch.abs(temperature - _METHANE_TC) / _METHANE_TC
    reduced_density_distance = (
        torch.abs(mass_density_g_cm3 - _CRITICAL_DENSITY_G_CM3) / _CRITICAL_DENSITY_G_CM3
    )
    damping = torch.exp(
        -18.66 * reduced_temperature_distance.square() - 4.25 * reduced_density_distance.pow(4)
    )
    return (
        prefactor
        * _BOLTZMANN_CONSTANT
        * temperature.square()
        / (6.0 * torch.pi * viscosity * correlation_length)
        * pressure_temperature.square()
        * torch.sqrt(compressibility)
        * damping
    )

methane_thermal_conductivity

methane_thermal_conductivity(
    temperature: Tensor,
    density_mol_l: Tensor,
    *,
    include_critical_enhancement: bool = True,
) -> Tensor

Evaluate the Hanley methane reference thermal conductivity.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Methane temperature in K.	required
density_mol_l	Tensor	Methane molar density in mol/L.	required
include_critical_enhancement	bool	Include the pure-fluid critical anomaly from Hanley Appendix C.	True

Returns:

Type	Description
Tensor	Thermal conductivity in W/(m K).

Notes

The background implements Pedersen et al. (2024), Eqs. 10.69-10.74, doi:10.1201/9780429457418-10. GT(9)=5.311764e-2 follows the primary Hanley table, doi:10.1016/0011-2275(75)90010-7. The critical term is a pure-methane correction and is excluded from the mixture model as directed by Pedersen section 10.2.

Source code in src/torch_flash/transport/thermal_conductivity.py

def methane_thermal_conductivity(
    temperature: Tensor,
    density_mol_l: Tensor,
    *,
    include_critical_enhancement: bool = True,
) -> Tensor:
    """Evaluate the Hanley methane reference thermal conductivity.

    Parameters
    ----------
    temperature
        Methane temperature in K.
    density_mol_l
        Methane molar density in mol/L.
    include_critical_enhancement
        Include the pure-fluid critical anomaly from Hanley Appendix C.

    Returns
    -------
    Tensor
        Thermal conductivity in W/(m K).

    Notes
    -----
    The background implements Pedersen et al. (2024), Eqs. 10.69-10.74,
    doi:10.1201/9780429457418-10. ``GT(9)=5.311764e-2`` follows the primary
    Hanley table, doi:10.1016/0011-2275(75)90010-7. The critical term is a
    pure-methane correction and is excluded from the mixture model as directed
    by Pedersen section 10.2.
    """
    temperature, density = torch.broadcast_tensors(temperature, density_mol_l)
    if bool(
        (
            (~torch.isfinite(temperature))
            | (~torch.isfinite(density))
            | (temperature <= 0.0)
            | (density < 0.0)
        ).any()
    ):
        raise InvalidStateError("methane conductivity requires positive T and non-negative density")
    background = _methane_conductivity_background(temperature, density)
    if not include_critical_enhancement or bool((density == 0.0).all()):
        return background
    enhancement = torch.where(
        density > 0.0,
        methane_critical_thermal_conductivity_enhancement(
            temperature,
            torch.clamp_min(density, torch.finfo(density.dtype).tiny),
        ),
        torch.zeros_like(density),
    )
    return background + enhancement

corresponding_states_viscosity

corresponding_states_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    components: ComponentSet,
    *,
    phase: Literal["liquid", "vapor"] = "vapor",
) -> Tensor

Evaluate the Pedersen corresponding-states mixture viscosity.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Positive temperature in K with arbitrary leading batch dimensions.	required
pressure	Tensor	Positive pressure in Pa, broadcastable with the state/composition batch.	required
composition	Tensor	Mole fractions with components on the final axis.	required
components	ComponentSet	Ordered critical-property and molar-mass data matching composition.	required
phase	Literal['liquid', 'vapor']	Methane-reference density root, "liquid" or "vapor".	'vapor'

Returns:

Type	Description
Tensor	Mixture dynamic viscosity in Pa s with the broadcast batch shape.

Raises:

Type	Description
ValueError	If the composition has no component axis or its component count differs from components.
InvalidStateError	If a required methane reference-state density cannot be solved.

Notes

This is a homogeneous-phase correlation and does not perform phase equilibrium. Supply the composition and phase root returned by the relevant state or flash calculation.

Source code in src/torch_flash/transport/viscosity.py

def corresponding_states_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    components: ComponentSet,
    *,
    phase: Literal["liquid", "vapor"] = "vapor",
) -> Tensor:
    """Evaluate the Pedersen corresponding-states mixture viscosity.

    Parameters
    ----------
    temperature
        Positive temperature in K with arbitrary leading batch dimensions.
    pressure
        Positive pressure in Pa, broadcastable with the state/composition
        batch.
    composition
        Mole fractions with components on the final axis.
    components
        Ordered critical-property and molar-mass data matching ``composition``.
    phase
        Methane-reference density root, ``"liquid"`` or ``"vapor"``.

    Returns
    -------
    Tensor
        Mixture dynamic viscosity in Pa s with the broadcast batch shape.

    Raises
    ------
    ValueError
        If the composition has no component axis or its component count
        differs from ``components``.
    InvalidStateError
        If a required methane reference-state density cannot be solved.

    Notes
    -----
    This is a homogeneous-phase correlation and does not perform phase
    equilibrium. Supply the composition and phase root returned by the
    relevant state or flash calculation.
    """
    if composition.ndim < 1:
        raise ValueError("composition must have a final component axis")
    if composition.shape[-1] != components.ncomponents:
        raise ValueError("composition and component set sizes must match")
    z = normalize_composition(composition)
    tc_mix, pc_mix = _mixture_pseudocritical(components, z)
    mapped_t = temperature * _METHANE_TC / tc_mix
    mapped_p = pressure * _METHANE_PC / pc_mix
    reference_density = methane_bwr_density(mapped_t, mapped_p, phase=phase)
    reference_mass_density = reference_density * _METHANE_MOLAR_MASS_G / 1000.0
    reduced_density = reference_mass_density / _METHANE_CRITICAL_MASS_DENSITY
    mixture_molar_mass = _pedersen_molecular_weight(components.molar_mass, z)
    alpha_mix = 1.0 + 7.378e-3 * reduced_density.pow(1.847) * mixture_molar_mass.pow(0.5173)
    alpha_reference = 1.0 + 7.378e-3 * reduced_density.pow(1.847) * temperature.new_tensor(
        _METHANE_MOLAR_MASS_G
    ).pow(0.5173)
    reference_pressure = pressure * _METHANE_PC * alpha_reference / (pc_mix * alpha_mix)
    reference_temperature = temperature * _METHANE_TC * alpha_reference / (tc_mix * alpha_mix)
    density = methane_bwr_density(reference_temperature, reference_pressure, phase=phase)
    reference_viscosity = methane_viscosity(reference_temperature, density)
    return (
        (tc_mix / _METHANE_TC).pow(-1.0 / 6.0)
        * (pc_mix / _METHANE_PC).pow(2.0 / 3.0)
        * (mixture_molar_mass / _METHANE_MOLAR_MASS_G).sqrt()
        * (alpha_mix / alpha_reference)
        * reference_viscosity
    )

friction_theory_viscosity

friction_theory_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    eos: CubicEOS,
    *,
    phase: PhaseKind = "stable",
    critical_viscosity: Tensor | None = None,
    critical_volume: Tensor | None = None,
) -> Tensor

Evaluate the one-parameter friction-theory cubic-EOS viscosity.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Positive state variables in K and Pa.	required
pressure	Tensor	Positive state variables in K and Pa.	required
composition	Tensor	Nonpolar-fluid mole fractions with components on the final axis.	required
eos	CubicEOS	Native SRK or PR76 :class:~torch_flash.eos.cubic.CubicEOS. Its mixing rule and binary interactions also define the pressure split.	required
phase	PhaseKind	Cubic root selection for the supplied homogeneous state.	'stable'
critical_viscosity	Tensor | None	Optional characteristic component viscosities in Pa s. When omitted, the n-alkane correlation of Eq. 21 is used.	None
critical_volume	Tensor | None	Optional component critical volumes in m3/mol for the Chung dilute-gas term. Defaults to the component database.	None

Returns:

Type	Description
Tensor	Dynamic viscosity in Pa s for each broadcast state.

Raises:

Type	Description
ValueError	If the EOS is not the published SRK or PR76 family, or component vectors have inconsistent shapes.
InvalidStateError	If a required volume or state is nonphysical or the result is nonpositive.

Notes

Implements Quiñones-Cisneros, Zéberg-Mikkelsen, and Stenby, "One parameter friction theory models for viscosity," Fluid Phase Equilibria 178 (2001), Eqs. 1-30 and Tables 1-2, doi:10.1016/S0378-3812(00)00474-X. The parameterization was developed for nonpolar fluids and n-alkane mixtures. PR78 and PRSV alpha functions are not silently identified with the paper's PR76 model.

Source code in src/torch_flash/transport/viscosity.py

def friction_theory_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    eos: CubicEOS,
    *,
    phase: PhaseKind = "stable",
    critical_viscosity: Tensor | None = None,
    critical_volume: Tensor | None = None,
) -> Tensor:
    """Evaluate the one-parameter friction-theory cubic-EOS viscosity.

    Parameters
    ----------
    temperature, pressure
        Positive state variables in K and Pa.
    composition
        Nonpolar-fluid mole fractions with components on the final axis.
    eos
        Native SRK or PR76 :class:`~torch_flash.eos.cubic.CubicEOS`. Its
        mixing rule and binary interactions also define the pressure split.
    phase
        Cubic root selection for the supplied homogeneous state.
    critical_viscosity
        Optional characteristic component viscosities in Pa s. When omitted,
        the n-alkane correlation of Eq. 21 is used.
    critical_volume
        Optional component critical volumes in m3/mol for the Chung dilute-gas
        term. Defaults to the component database.

    Returns
    -------
    Tensor
        Dynamic viscosity in Pa s for each broadcast state.

    Raises
    ------
    ValueError
        If the EOS is not the published SRK or PR76 family, or component
        vectors have inconsistent shapes.
    InvalidStateError
        If a required volume or state is nonphysical or the result is
        nonpositive.

    Notes
    -----
    Implements Quiñones-Cisneros, Zéberg-Mikkelsen, and Stenby,
    "One parameter friction theory models for viscosity," *Fluid Phase
    Equilibria* 178 (2001), Eqs. 1-30 and Tables 1-2,
    doi:10.1016/S0378-3812(00)00474-X. The parameterization was developed for
    nonpolar fluids and n-alkane mixtures. PR78 and PRSV alpha functions are
    not silently identified with the paper's PR76 model.
    """
    if eos.constants.alpha_kind == "srk":
        family: Literal["SRK", "PR"] = "SRK"
    elif eos.constants.alpha_kind == "pr76":
        family = "PR"
    else:
        raise ValueError("friction theory supports only the published SRK and PR76 families")
    if composition.ndim < 1 or composition.shape[-1] != eos.ncomponents:
        raise ValueError("composition and cubic EOS component sizes must match")
    if bool(
        (
            (~torch.isfinite(temperature))
            | (~torch.isfinite(pressure))
            | (temperature <= 0.0)
            | (pressure <= 0.0)
        ).any()
    ):
        raise InvalidStateError(
            "friction-theory temperature and pressure must be finite and positive"
        )
    z = normalize_composition(composition)
    volumes = eos.critical_volume if critical_volume is None else critical_volume
    if volumes is None or volumes.shape != eos.critical_temperature.shape:
        raise ValueError("friction theory requires one critical volume per component")
    if bool(((~torch.isfinite(volumes)) | (volumes <= 0.0)).any()):
        raise ValueError("critical volumes must be finite and positive")
    if critical_viscosity is None:
        pressure_bar = eos.critical_pressure / 1.0e5
        molecular_weight_g = 1000.0 * eos.molar_mass
        eta_c = 1.0e-7 * 0.597556 * pressure_bar * molecular_weight_g.pow(0.601652)
    else:
        eta_c = critical_viscosity.to(dtype=temperature.dtype, device=temperature.device)
        if eta_c.shape != eos.critical_temperature.shape:
            raise ValueError("critical_viscosity must have one value per component")
        if bool(((~torch.isfinite(eta_c)) | (eta_c <= 0.0)).any()):
            raise ValueError("critical viscosities must be finite and positive")
    kr_hat, ka_hat, krr_hat = _friction_reduced_coefficients(
        temperature,
        eos,
        family,
    )
    molecular_weight_g = 1000.0 * eos.molar_mass
    epsilon = 0.30
    denominator = torch.sum(z / molecular_weight_g.pow(epsilon), dim=-1, keepdim=True)
    weights = z / (molecular_weight_g.pow(epsilon) * denominator)
    kr = torch.sum(weights * eta_c * kr_hat / eos.critical_pressure, dim=-1)
    ka = torch.sum(weights * eta_c * ka_hat / eos.critical_pressure, dim=-1)
    krr = torch.sum(weights * eta_c * krr_hat / eos.critical_pressure.square(), dim=-1)
    dilute_components = _chung_dilute_viscosity(temperature, eos, volumes)
    dilute = torch.exp(torch.sum(z * torch.log(dilute_components), dim=-1))
    physical_volume = eos.molar_volume(temperature, pressure, z, phase)
    translation = torch.sum(
        z * eos.component_volume_translation(temperature),
        dim=-1,
    )
    volume = physical_volume - translation
    attraction, covolume = eos.mixture_parameters(temperature, z)
    repulsive_pressure = R * temperature / (volume - covolume)
    attractive_pressure = -attraction / (
        (volume + eos.constants.delta1 * covolume) * (volume + eos.constants.delta2 * covolume)
    )
    viscosity = (
        dilute
        + kr * repulsive_pressure
        + ka * attractive_pressure
        + krr * repulsive_pressure.square()
    )
    if bool((~torch.isfinite(viscosity) | (viscosity <= 0.0)).any()):
        raise InvalidStateError("friction theory produced a non-positive viscosity")
    return viscosity

heavy_oil_corresponding_states_viscosity

heavy_oil_corresponding_states_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    components: ComponentSet,
    *,
    phase: Literal["liquid", "vapor"] = "liquid",
    third_csp: Tensor | float = 1.0,
    fourth_csp: Tensor | float = 1.0,
) -> Tensor

Evaluate the blended corresponding-states heavy-oil viscosity model.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Positive SI state variables in K and Pa.	required
pressure	Tensor	Positive SI state variables in K and Pa.	required
composition	Tensor	Mole fractions with components on the final axis.	required
components	ComponentSet	Critical properties and molar masses in the same order.	required
phase	Literal['liquid', 'vapor']	Methane-reference BWR density branch.	'liquid'
third_csp	Tensor | float	Dimensionless Lindeloff tuning factors.	1.0
fourth_csp	Tensor | float	Dimensionless Lindeloff tuning factors.	1.0

Returns:

Type	Description
Tensor	Dynamic viscosity in Pa s over the broadcast state/composition batch.

Notes

Pedersen et al. (2024), section 10.1.2, retains the conventional corresponding-states result above a methane reference temperature of 75 K, uses the Lindeloff branch below 65 K, and linearly blends the two between 65 and 75 K. Calibration factors are explicit tensors so fitting preserves PyTorch gradients.

Source code in src/torch_flash/transport/viscosity.py

def heavy_oil_corresponding_states_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    composition: Tensor,
    components: ComponentSet,
    *,
    phase: Literal["liquid", "vapor"] = "liquid",
    third_csp: Tensor | float = 1.0,
    fourth_csp: Tensor | float = 1.0,
) -> Tensor:
    """Evaluate the blended corresponding-states heavy-oil viscosity model.

    Parameters
    ----------
    temperature, pressure
        Positive SI state variables in K and Pa.
    composition
        Mole fractions with components on the final axis.
    components
        Critical properties and molar masses in the same order.
    phase
        Methane-reference BWR density branch.
    third_csp, fourth_csp
        Dimensionless Lindeloff tuning factors.

    Returns
    -------
    Tensor
        Dynamic viscosity in Pa s over the broadcast state/composition batch.

    Notes
    -----
    Pedersen et al. (2024), section 10.1.2, retains the conventional
    corresponding-states result above a methane reference temperature of
    75 K, uses the Lindeloff branch below 65 K, and linearly blends the two
    between 65 and 75 K. Calibration factors are explicit tensors so fitting
    preserves PyTorch gradients.
    """
    if composition.ndim < 1 or composition.shape[-1] != components.ncomponents:
        raise ValueError("composition and component set sizes must match")
    z = normalize_composition(composition)
    tc_mix, pc_mix = _mixture_pseudocritical(components, z)
    mapped_temperature = temperature * _METHANE_TC / tc_mix
    mapped_pressure = pressure * _METHANE_PC / pc_mix
    initial_density = methane_bwr_density(mapped_temperature, mapped_pressure, phase=phase)
    reduced_density = (
        initial_density * _METHANE_MOLAR_MASS_G / 1000.0 / _METHANE_CRITICAL_MASS_DENSITY
    )
    molar_mass_g = 1000.0 * components.molar_mass
    mn_g = torch.sum(z * molar_mass_g, dim=-1)
    mw_g = torch.sum(z * molar_mass_g.square(), dim=-1) / mn_g
    mixture_mass = 1.304e-4 * (mw_g.pow(2.303) - mn_g.pow(2.303)) + mn_g
    alpha_mix = 1.0 + 7.378e-3 * reduced_density.pow(1.847) * mixture_mass.pow(0.5173)
    alpha_reference = 1.0 + 7.378e-3 * reduced_density.pow(1.847) * temperature.new_tensor(
        _METHANE_MOLAR_MASS_G
    ).pow(0.5173)
    reference_temperature = temperature * _METHANE_TC * alpha_reference / (tc_mix * alpha_mix)
    conventional = corresponding_states_viscosity(
        temperature,
        pressure,
        z,
        components,
        phase=phase,
    )
    heavy = stabilized_heavy_oil_viscosity(
        temperature,
        pressure,
        mn_g / 1000.0,
        mw_g / 1000.0,
        third_csp=third_csp,
        fourth_csp=fourth_csp,
    )
    heavy_weight = torch.clamp((75.0 - reference_temperature) / 10.0, 0.0, 1.0)
    return torch.lerp(conventional, heavy, heavy_weight)

kinematic_viscosity

kinematic_viscosity(
    dynamic_viscosity: Tensor, mass_density: Tensor
) -> Tensor

Return kinematic viscosity from dynamic viscosity and mass density.

Parameters:

Name	Type	Description	Default
dynamic_viscosity	Tensor	Dynamic viscosity in Pa s.	required
mass_density	Tensor	Homogeneous-phase mass density in kg/m3, broadcastable with dynamic_viscosity.	required

Returns:

Type	Description
Tensor	Kinematic viscosity in m2/s with the broadcast input shape.

Raises:

Type	Description
InvalidStateError	If viscosity is negative or either input is nonfinite or density is nonpositive.

Notes

This is the SI definition in Pedersen et al. (2024), section 10.1, rather than an equilibrium or empirical transport model.

Source code in src/torch_flash/transport/viscosity.py

def kinematic_viscosity(dynamic_viscosity: Tensor, mass_density: Tensor) -> Tensor:
    """Return kinematic viscosity from dynamic viscosity and mass density.

    Parameters
    ----------
    dynamic_viscosity
        Dynamic viscosity in Pa s.
    mass_density
        Homogeneous-phase mass density in kg/m3, broadcastable with
        ``dynamic_viscosity``.

    Returns
    -------
    Tensor
        Kinematic viscosity in m2/s with the broadcast input shape.

    Raises
    ------
    InvalidStateError
        If viscosity is negative or either input is nonfinite or density is
        nonpositive.

    Notes
    -----
    This is the SI definition in Pedersen et al. (2024), section 10.1, rather
    than an equilibrium or empirical transport model.
    """
    viscosity, density = torch.broadcast_tensors(dynamic_viscosity, mass_density)
    invalid = (
        (~torch.isfinite(viscosity))
        | (~torch.isfinite(density))
        | (viscosity < 0.0)
        | (density <= 0.0)
    )
    if bool(invalid.any()):
        raise InvalidStateError("viscosity must be non-negative and density finite and positive")
    result: Tensor = viscosity / density
    return result

lbc_pseudocomponent_critical_volume

lbc_pseudocomponent_critical_volume(
    molar_mass: Tensor, standard_liquid_density: Tensor
) -> Tensor

Estimate a C7+ critical molar volume from Pedersen Eq. 10.41.

Parameters:

Name	Type	Description	Default
molar_mass	Tensor	Pseudo-component molar mass in kg/mol.	required
standard_liquid_density	Tensor	Standard liquid mass density in kg/m3.	required

Returns:

Type	Description
Tensor	Estimated critical molar volume in m3/mol.

Raises:

Type	Description
ValueError	If either input is nonfinite or nonpositive.

Notes

Parameters are SI: molar mass in kg/mol and liquid density in kg/m3. The returned critical molar volume is in m3/mol. The underlying empirical equation uses g/mol, g/cm3, and ft3/lbmol.

Source code in src/torch_flash/transport/viscosity.py

def lbc_pseudocomponent_critical_volume(
    molar_mass: Tensor,
    standard_liquid_density: Tensor,
) -> Tensor:
    """Estimate a C7+ critical molar volume from Pedersen Eq. 10.41.

    Parameters
    ----------
    molar_mass
        Pseudo-component molar mass in kg/mol.
    standard_liquid_density
        Standard liquid mass density in kg/m3.

    Returns
    -------
    Tensor
        Estimated critical molar volume in m3/mol.

    Raises
    ------
    ValueError
        If either input is nonfinite or nonpositive.

    Notes
    -----
    Parameters are SI: molar mass in kg/mol and liquid density in kg/m3.
    The returned critical molar volume is in m3/mol. The underlying empirical
    equation uses g/mol, g/cm3, and ft3/lbmol.
    """
    if bool(
        (
            (~torch.isfinite(molar_mass))
            | (~torch.isfinite(standard_liquid_density))
            | (molar_mass <= 0.0)
            | (standard_liquid_density <= 0.0)
        ).any()
    ):
        raise ValueError("molar mass and standard liquid density must be finite and positive")
    molar_mass_g = 1000.0 * molar_mass
    density_g_cm3 = standard_liquid_density / 1000.0
    volume_ft3_lbmol = (
        21.573
        + 0.015122 * molar_mass_g
        - 27.656 * density_g_cm3
        + 0.070615 * molar_mass_g * density_g_cm3
    )
    return volume_ft3_lbmol * _FT3_PER_LBMOL_TO_M3_PER_MOL

lbc_viscosity

lbc_viscosity(
    temperature: Tensor,
    molar_density: Tensor,
    composition: Tensor,
    components: ComponentSet,
    *,
    critical_volume: Tensor | None = None,
    coefficients: Tensor | None = None,
) -> Tensor

Return Lohrenz-Bray-Clark mixture viscosity in Pa s.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperature in K.	required
molar_density	Tensor	Homogeneous-phase molar density in mol/m3, normally obtained from the same EoS used by the flash calculation.	required
composition	Tensor	Mole fractions with components on the final axis.	required
components	ComponentSet	Critical temperatures, pressures, molar masses, and volumes.	required
critical_volume	Tensor | None	Optional m3/mol values overriding the component data. This is the principal LBC tuning parameter for petroleum pseudocomponents.	None
coefficients	Tensor | None	Optional trainable or fitted a1...a5 tensor replacing the original LBC constants.	None

Returns:

Type	Description
Tensor	Mixture dynamic viscosity in Pa s.

Raises:

Type	Description
ValueError	If component shapes, critical volumes, or coefficient tensors are inconsistent.
InvalidStateError	If temperature/density or a required reducing property is outside its physical domain.

Notes

The Stiel-Thodos reducing parameters use K, atm, and g/mol, as required by the published numerical constants. LBC is empirical and often needs pseudocomponent critical-volume or coefficient tuning for heavy oils.

Source code in src/torch_flash/transport/viscosity.py

def lbc_viscosity(
    temperature: Tensor,
    molar_density: Tensor,
    composition: Tensor,
    components: ComponentSet,
    *,
    critical_volume: Tensor | None = None,
    coefficients: Tensor | None = None,
) -> Tensor:
    """Return Lohrenz-Bray-Clark mixture viscosity in Pa s.

    Parameters
    ----------
    temperature
        Temperature in K.
    molar_density
        Homogeneous-phase molar density in mol/m3, normally obtained from the
        same EoS used by the flash calculation.
    composition
        Mole fractions with components on the final axis.
    components
        Critical temperatures, pressures, molar masses, and volumes.
    critical_volume
        Optional m3/mol values overriding the component data. This is the
        principal LBC tuning parameter for petroleum pseudocomponents.
    coefficients
        Optional trainable or fitted ``a1...a5`` tensor replacing the original
        LBC constants.

    Returns
    -------
    Tensor
        Mixture dynamic viscosity in Pa s.

    Raises
    ------
    ValueError
        If component shapes, critical volumes, or coefficient tensors are
        inconsistent.
    InvalidStateError
        If temperature/density or a required reducing property is outside its
        physical domain.

    Notes
    -----
    The Stiel-Thodos reducing parameters use K, atm, and g/mol, as required by
    the published numerical constants. LBC is empirical and often needs
    pseudocomponent critical-volume or coefficient tuning for heavy oils.
    """
    if composition.ndim < 1:
        raise ValueError("composition must have a final component axis")
    if composition.shape[-1] != components.ncomponents:
        raise ValueError("composition and component set sizes must match")
    if bool(
        (
            (~torch.isfinite(temperature))
            | (~torch.isfinite(molar_density))
            | (temperature <= 0.0)
            | (molar_density < 0.0)
        ).any()
    ):
        raise InvalidStateError("temperature must be positive and molar density non-negative")

    volumes = components.critical_volume if critical_volume is None else critical_volume
    if volumes is None:
        raise ValueError("critical molar volumes are required by the LBC correlation")
    if volumes.shape != components.critical_temperature.shape:
        raise ValueError("critical_volume must have one value per component")
    if bool(((~torch.isfinite(volumes)) | (volumes <= 0.0)).any()):
        raise ValueError("critical molar volumes must be finite and positive")

    parameters = (
        temperature.new_tensor(_LBC)
        if coefficients is None
        else coefficients.to(dtype=temperature.dtype, device=temperature.device)
    )
    if parameters.shape != (5,) or bool((~torch.isfinite(parameters)).any()):
        raise ValueError("LBC coefficients must be a finite five-element tensor")

    z = normalize_composition(composition)
    tc = components.critical_temperature
    pc_atm = components.critical_pressure / 101_325.0
    molar_mass_g = 1000.0 * components.molar_mass
    xi_components = tc.pow(1.0 / 6.0) / (molar_mass_g.sqrt() * pc_atm.pow(2.0 / 3.0))
    reduced_temperature = temperature[..., None] / tc
    dilute_components_cp = torch.where(
        reduced_temperature <= 1.5,
        34.0e-5 * reduced_temperature.pow(0.94) / xi_components,
        17.78e-5 * (4.58 * reduced_temperature - 1.67).pow(5.0 / 8.0) / xi_components,
    )
    square_root_mass = molar_mass_g.sqrt()
    dilute_mixture_cp = torch.sum(
        z * dilute_components_cp * square_root_mass,
        dim=-1,
    ) / torch.sum(z * square_root_mass, dim=-1)

    mixture_xi = torch.sum(z * tc, dim=-1).pow(1.0 / 6.0) / (
        torch.sum(z * molar_mass_g, dim=-1).sqrt() * torch.sum(z * pc_atm, dim=-1).pow(2.0 / 3.0)
    )
    reduced_density = molar_density * torch.sum(z * volumes, dim=-1)
    polynomial = parameters[0] + reduced_density * (
        parameters[1]
        + reduced_density
        * (parameters[2] + reduced_density * (parameters[3] + reduced_density * parameters[4]))
    )
    dense_increment_cp = (polynomial.pow(4) - 1.0e-4) / mixture_xi
    viscosity_cp = dilute_mixture_cp + dense_increment_cp
    if bool((~torch.isfinite(viscosity_cp) | (viscosity_cp <= 0.0)).any()):
        raise InvalidStateError("LBC correlation produced a non-positive viscosity")
    return viscosity_cp * 1.0e-3

lee_gas_viscosity

lee_gas_viscosity(
    temperature: Tensor,
    mass_density: Tensor,
    molar_mass: Tensor,
) -> Tensor

Evaluate the Lee-Gonzalez-Eakin natural-gas viscosity correlation.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Gas temperature in K.	required
mass_density	Tensor	Gas mass density in kg/m3.	required
molar_mass	Tensor	Mixture molar mass in kg/mol.	required

Returns:

Type	Description
Tensor	Dynamic viscosity in Pa s with the broadcast input shape.

Raises:

Type	Description
InvalidStateError	If temperature or molar mass is nonpositive, density is negative, or any input or output is nonfinite.

Notes

Implements Lee, Gonzalez, and Eakin, "The Viscosity of Natural Gases," J. Petroleum Technology 18 (1966), 997-1000, doi:10.2118/1340-PA, as reproduced in Pedersen et al. (2024), Eqs. 10.46-10.49. The published equation uses degrees Rankine, lb/ft3, g/mol, and cP; this API converts SI inputs and returns SI viscosity. It is a gas correlation and is not valid for a condensed liquid phase.

Source code in src/torch_flash/transport/viscosity.py

def lee_gas_viscosity(
    temperature: Tensor,
    mass_density: Tensor,
    molar_mass: Tensor,
) -> Tensor:
    """Evaluate the Lee-Gonzalez-Eakin natural-gas viscosity correlation.

    Parameters
    ----------
    temperature
        Gas temperature in K.
    mass_density
        Gas mass density in kg/m3.
    molar_mass
        Mixture molar mass in kg/mol.

    Returns
    -------
    Tensor
        Dynamic viscosity in Pa s with the broadcast input shape.

    Raises
    ------
    InvalidStateError
        If temperature or molar mass is nonpositive, density is negative, or
        any input or output is nonfinite.

    Notes
    -----
    Implements Lee, Gonzalez, and Eakin, "The Viscosity of Natural Gases,"
    *J. Petroleum Technology* 18 (1966), 997-1000,
    doi:10.2118/1340-PA, as reproduced in Pedersen et al. (2024),
    Eqs. 10.46-10.49. The published equation uses degrees Rankine, lb/ft3,
    g/mol, and cP; this API converts SI inputs and returns SI viscosity. It is
    a gas correlation and is not valid for a condensed liquid phase.
    """
    temperature, density, molecular_weight = torch.broadcast_tensors(
        temperature,
        mass_density,
        molar_mass,
    )
    invalid = (
        (~torch.isfinite(temperature))
        | (~torch.isfinite(density))
        | (~torch.isfinite(molecular_weight))
        | (temperature <= 0.0)
        | (density < 0.0)
        | (molecular_weight <= 0.0)
    )
    if bool(invalid.any()):
        raise InvalidStateError(
            "temperature and molar mass must be positive and gas density non-negative"
        )
    temperature_rankine = 1.8 * temperature
    molecular_weight_g = 1000.0 * molecular_weight
    density_lb_ft3 = density / _LB_FT3_TO_KG_M3
    x = 3.5 + 986.0 / temperature_rankine + 0.01 * molecular_weight_g
    y = 2.4 - 0.2 * x
    k = (
        (9.4 + 0.02 * molecular_weight_g)
        * temperature_rankine.pow(1.5)
        / (209.0 + 19.0 * molecular_weight_g + temperature_rankine)
    )
    viscosity_cp = 1.0e-4 * k * torch.exp(x * (density_lb_ft3 / 62.4).pow(y))
    if bool((~torch.isfinite(viscosity_cp) | (viscosity_cp <= 0.0)).any()):
        raise InvalidStateError("Lee gas correlation produced a non-positive viscosity")
    result: Tensor = 1.0e-3 * viscosity_cp
    return result

methane_viscosity

methane_viscosity(
    temperature: Tensor, density_mol_l: Tensor
) -> Tensor

Return methane dynamic viscosity in Pa s from Eq. 10.6.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperature in K.	required
density_mol_l	Tensor	Methane molar density in mol/L.	required

Returns:

Type	Description
Tensor	Dynamic viscosity in Pa s.

Notes

Although the 2024 textbook prose labels the correlation density as mol/L, the published Hanley coefficients require mass density in kg/L (numerically equal to g/cm3). Using molar density inside the fractional powers produces viscosities two orders of magnitude too large at reservoir pressures.

Source code in src/torch_flash/transport/viscosity.py

def methane_viscosity(
    temperature: Tensor,
    density_mol_l: Tensor,
) -> Tensor:
    """Return methane dynamic viscosity in Pa s from Eq. 10.6.

    Parameters
    ----------
    temperature
        Temperature in K.
    density_mol_l
        Methane molar density in mol/L.

    Returns
    -------
    Tensor
        Dynamic viscosity in Pa s.

    Notes
    -----
    Although the 2024 textbook prose labels the correlation density as mol/L,
    the published Hanley coefficients require mass density in kg/L (numerically
    equal to g/cm3). Using molar density inside the fractional powers produces
    viscosities two orders of magnitude too large at reservoir pressures.
    """
    gv = torch.tensor(_GV, dtype=temperature.dtype, device=temperature.device)
    powers = torch.arange(-3, 6, dtype=temperature.dtype, device=temperature.device) / 3.0
    dilute = torch.sum(gv * temperature[..., None].pow(powers), dim=-1)
    eta1 = 1.696985927 - 0.133372346 * (1.4 - torch.log(temperature / 168.0)).square()
    mass_density = density_mol_l * _METHANE_MOLAR_MASS_G / 1000.0
    theta = (mass_density - _METHANE_CRITICAL_MASS_DENSITY) / _METHANE_CRITICAL_MASS_DENSITY
    j = torch.tensor(_J, dtype=temperature.dtype, device=temperature.device)
    dense = torch.exp(j[0] + j[3] / temperature) * (
        torch.exp(
            mass_density.pow(0.1) * (j[1] + j[2] / temperature.pow(1.5))
            + theta
            * mass_density.sqrt()
            * (j[4] + j[5] / temperature + j[6] / temperature.square())
        )
        - 1.0
    )
    low_values = cast(Mapping[str, object], _TRANSPORT_PARAMETERS["methane"])[
        "viscosity_low_temperature"
    ]
    low = torch.as_tensor(low_values, dtype=temperature.dtype, device=temperature.device)
    low_dense = torch.exp(low[0] + low[3] / temperature) * (
        torch.exp(
            mass_density.pow(0.1) * (low[1] + low[2] / temperature.pow(1.5))
            + mass_density.sqrt() * (low[4] + low[5] / temperature + low[6] / temperature.square())
        )
        - 1.0
    )
    transition = 0.5 * (torch.tanh(temperature - 91.0) + 1.0)
    dense = transition * dense + (1.0 - transition) * low_dense
    viscosity_cp = 1.0e-4 * (dilute + eta1 * mass_density + dense)
    return viscosity_cp * 1.0e-3

stabilized_heavy_oil_viscosity

stabilized_heavy_oil_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    number_average_molar_mass: Tensor,
    weight_average_molar_mass: Tensor,
    *,
    third_csp: Tensor | float = 1.0,
    fourth_csp: Tensor | float = 1.0,
) -> Tensor

Evaluate the Lindeloff stabilized/live-heavy-oil viscosity branch.

Parameters:

Name	Type	Description	Default
temperature	Tensor	Temperature in K.	required
pressure	Tensor	Pressure in Pa.	required
number_average_molar_mass	Tensor	Mixture molecular-weight averages in kg/mol.	required
weight_average_molar_mass	Tensor	Mixture molecular-weight averages in kg/mol.	required
third_csp	Tensor | float	Dimensionless fitted factors in Pedersen et al. (2024), Eqs. 10.34-10.36. Defaults reproduce the published predictive branch.	1.0
fourth_csp	Tensor | float	Dimensionless fitted factors in Pedersen et al. (2024), Eqs. 10.34-10.36. Defaults reproduce the published predictive branch.	1.0

Returns:

Type	Description
Tensor	Dynamic viscosity in Pa s.

Raises:

Type	Description
InvalidStateError	If the state or molecular-weight averages are outside their physical domains, or the correlation returns a nonfinite result.

Notes

Implements Lindeloff et al., "The corresponding states viscosity model applied to heavy oil systems," J. Can. Petroleum Technology 43 (2004), 47-53, through Pedersen et al. (2024), Eqs. 10.33-10.37. This function is the empirical heavy-oil branch only; use :func:heavy_oil_corresponding_states_viscosity for the published low-reference-temperature blending protocol.

Source code in src/torch_flash/transport/viscosity.py

def stabilized_heavy_oil_viscosity(
    temperature: Tensor,
    pressure: Tensor,
    number_average_molar_mass: Tensor,
    weight_average_molar_mass: Tensor,
    *,
    third_csp: Tensor | float = 1.0,
    fourth_csp: Tensor | float = 1.0,
) -> Tensor:
    """Evaluate the Lindeloff stabilized/live-heavy-oil viscosity branch.

    Parameters
    ----------
    temperature
        Temperature in K.
    pressure
        Pressure in Pa.
    number_average_molar_mass, weight_average_molar_mass
        Mixture molecular-weight averages in kg/mol.
    third_csp, fourth_csp
        Dimensionless fitted factors in Pedersen et al. (2024),
        Eqs. 10.34-10.36. Defaults reproduce the published predictive branch.

    Returns
    -------
    Tensor
        Dynamic viscosity in Pa s.

    Raises
    ------
    InvalidStateError
        If the state or molecular-weight averages are outside their physical
        domains, or the correlation returns a nonfinite result.

    Notes
    -----
    Implements Lindeloff et al., "The corresponding states viscosity model
    applied to heavy oil systems," *J. Can. Petroleum Technology* 43 (2004),
    47-53, through Pedersen et al. (2024), Eqs. 10.33-10.37. This function is
    the empirical heavy-oil branch only; use
    :func:`heavy_oil_corresponding_states_viscosity` for the published
    low-reference-temperature blending protocol.
    """
    third = torch.as_tensor(third_csp, dtype=temperature.dtype, device=temperature.device)
    fourth = torch.as_tensor(fourth_csp, dtype=temperature.dtype, device=temperature.device)
    temperature, pressure, mn, mw, third, fourth = torch.broadcast_tensors(
        temperature,
        pressure,
        number_average_molar_mass,
        weight_average_molar_mass,
        third,
        fourth,
    )
    invalid = (
        (~torch.isfinite(temperature))
        | (~torch.isfinite(pressure))
        | (~torch.isfinite(mn))
        | (~torch.isfinite(mw))
        | (~torch.isfinite(third))
        | (~torch.isfinite(fourth))
        | (temperature <= 0.0)
        | (pressure <= 0.0)
        | (mn <= 0.0)
        | (mw < mn)
        | (third <= 0.0)
        | (fourth <= 0.0)
    )
    if bool(invalid.any()):
        raise InvalidStateError(
            "heavy-oil state, molecular-weight averages, and CSP factors must be physical"
        )
    mn_g = 1000.0 * mn
    mw_g = 1000.0 * mw
    visfac3 = 0.2252 * temperature / mn_g + 0.9738
    visfac4 = 0.5354 * visfac3 - 0.1170
    ratio = mw_g / mn_g
    low_base = 1.5 / (visfac3 * third)
    high_base = mw_g / (visfac3 * third * mn_g)
    representative_mass = mn_g * torch.where(
        ratio <= 1.5,
        low_base.pow(visfac4 * fourth),
        high_base.pow(visfac4 * fourth),
    )
    mass_sign = torch.where(temperature > 564.49, 1.0, -1.0)
    log10_viscosity_cp = (
        -0.07995
        + mass_sign * 0.01101 * representative_mass
        - 371.8 / temperature
        + 6.215 * representative_mass / temperature
    )
    viscosity_atmospheric_cp = torch.pow(temperature.new_tensor(10.0), log10_viscosity_cp)
    pressure_atm = pressure / 101_325.0
    pressure_factor = torch.exp(0.00384 * (pressure_atm.pow(0.8226) - 1.0) / 0.8226)
    viscosity = 1.0e-3 * viscosity_atmospheric_cp * pressure_factor
    if bool((~torch.isfinite(viscosity) | (viscosity <= 0.0)).any()):
        raise InvalidStateError("heavy-oil correlation produced a non-positive viscosity")
    return viscosity