Source code for chempy.debye_huckel

# -*- coding: utf-8 -*-
"""
This modules collects expressions for the Debye-Hückel constants.

"""
from __future__ import (absolute_import, division, print_function)

from .core import _ActivityProductBase, ionic_strength
from ._util import get_backend


[docs]def A(eps_r, T, rho, b0=1, constants=None, units=None, backend=None): """ Debye Huckel constant A Parameters ---------- eps_r: float relative permittivity T: float with unit Temperature (default: assume Kelvin) rho: float with unit density (default: assume kg/m**3) b0: float with unit reference molality (default: assume mol/kg) units: object (optional, default: None) attributes accessed: meter, Kelvin and mol constants: object (optional, default: None) if None: T assumed to be in Kelvin and b0 = 1 mol/kg else: see source code for what attributes are used. Tip: pass quantities.constants Notes ----- Remember to divide by ln(10) if you want to use the constant with log10 based expression. References ---------- Atkins, De Paula, Physical Chemistry, 8th edition """ be = get_backend(backend) one = be.pi**0 if constants is None: combined = 132871.85866393594 if units is not None: m = units.meter K = units.Kelvin mol = units.mol combined *= (m*K)**(3*one/2) / mol**(one/2) return combined*(rho * b0 * T**-3 * eps_r**-3)**0.5 F = constants.Faraday_constant NA = constants.Avogadro_constant eps0 = constants.vacuum_permittivity kB = constants.Boltzmann_constant pi = constants.pi A = F**3/(4*pi*NA)*(rho*b0/(2*(eps0*eps_r*kB*NA*T)**3))**(one/2) return A
[docs]def B(eps_r, T, rho, b0=1, constants=None, units=None, backend=None): """ Extended Debye-Huckel parameter B Parameters ---------- eps_r: float relative permittivity T: float with unit temperature rho: float with unit density b0: float with unit reference molality units: object (optional, default: None) attributes accessed: meter, Kelvin and mol constants: object (optional, default: None) if None: T assumed to be in Kelvin, rho in kg/m**3 and b0 = 1 mol/kg else: attributes accessed: molar_gas_constant, Faraday_constant Tip: pass quantities.constants Returns ------- Debye Huckel B constant (default in m**-1) """ be = get_backend(backend) one = be.pi**0 if constants is None: combined = 15903203868.740343 if units is not None: m = units.meter K = units.Kelvin mol = units.mol combined *= (m*K/mol)**(one/2) return combined*(rho*b0/(T*eps_r))**0.5 F = constants.Faraday_constant eps0 = constants.vacuum_permittivity R = constants.molar_gas_constant B = F*(2*rho*b0/(eps_r*eps0*R*T))**(one/2) return B
[docs]def limiting_log_gamma(I, z, A, I0=1, backend=None): """ Debye-Hyckel limiting formula """ be = get_backend(backend) one = be.pi**0 return -A*z**2*(I/I0)**(one/2)
[docs]def extended_log_gamma(I, z, a, A, B, C=0, I0=1, backend=None): """ Debye-Huckel extended formula """ be = get_backend(backend) one = be.pi**0 I_I0 = I/I0 sqrt_I_I0 = (I_I0)**(one/2) return -A*z**2 * sqrt_I_I0/(1 + B*a*sqrt_I_I0) + C*I_I0
[docs]def davies_log_gamma(I, z, A, C=-0.3, I0=1, backend=None): """ Davies formula """ be = get_backend(backend) one = be.pi**0 I_I0 = I/I0 sqrt_I_I0 = (I_I0)**(one/2) return -A * z**2 * (sqrt_I_I0/(1 + sqrt_I_I0) + C*I_I0)
[docs]def limiting_activity_product(I, stoich, z, T, eps_r, rho, backend=None): """ Product of activity coefficients based on DH limiting law. """ be = get_backend(backend) Aval = A(eps_r, T, rho) tot = 0 for idx, nr in enumerate(stoich): tot += nr*limiting_log_gamma(I, z[idx], Aval) return be.exp(tot)
[docs]def extended_activity_product(I, stoich, z, a, T, eps_r, rho, C=0, backend=None): be = get_backend(backend) Aval = A(eps_r, T, rho) Bval = B(eps_r, T, rho) tot = 0 for idx, nr in enumerate(stoich): tot += nr*extended_log_gamma(I, z[idx], a[idx], Aval, Bval, C) return be.exp(tot)
[docs]def davies_activity_product(I, stoich, z, a, T, eps_r, rho, C=-0.3, backend=None): be = get_backend(backend) Aval = A(eps_r, T, rho) tot = 0 for idx, nr in enumerate(stoich): tot += nr*davies_log_gamma(I, z[idx], Aval, C) return be.exp(tot)
[docs]class LimitingDebyeHuckelActivityProduct(_ActivityProductBase): def __call__(self, c): z = self.args[0] I = ionic_strength(c, z) return limiting_activity_product(I, self.stoich, *self.args)
[docs]class ExtendedDebyeHuckelActivityProduct(_ActivityProductBase): def __call__(self, c): z = self.args[0] I = ionic_strength(c, z) return extended_activity_product(I, self.stoich, *self.args)