Add porous media calculations and example

Author: Robert DeJaco

README.md

@@ -1,5 +1,5 @@
 # diffusivity_estimation
-Estimation of Diffusion Coefficients
+Estimation of Diffusion Coefficients for gases and porous media
 
 Diffusion coefficients are often estimated in scientific and engineering.
 There can be many options and parameters involved in these estimations.
@@ -8,3 +8,8 @@ Having an open-source software repository for estimation of diffusion coefficien
 streamline reproduction and sensitivity analysis.
 
 So far, this repository focuses on estimation of diffusivities in porous media.
+
+
+Binary Diffusivities approximate range (Deen p. 11)
+===================================================
+* Gases from 1E-5 to 1E-4

examples/methane_hydrogen_porous_media/input.csv

@@ -0,0 +1,17 @@
+void_fraction, ,0.35
+d_pore, ,5e-10
+molecular_weight, ,None
+tortuosity, ,2.0
+R, ,8.314
+molecular_weight_i, , 
+ ,CH4,16.043
+ ,H2,2.016
+sigma_i, , 
+ ,CH4,3.758
+ ,H2,2.827
+epsilon_i, , 
+ ,CH4,148.6
+ ,H2,59.7
+component_i, , 
+ ,CH4, 
+ ,H2, 

examples/methane_hydrogen_porous_media/output.csv

@@ -0,0 +1,6 @@
+CH4,H2,effective_macropore_i [m^2/s],1.015327e-10
+CH4,H2,knudsen_i [m^2/s],5.803363e-10
+CH4,H2,molecular_ij [m^2/s],2.220814e-11
+H2,CH4,effective_macropore_i [m^2/s],2.862856e-10
+H2,CH4,knudsen_i [m^2/s],1.637108e-09
+H2,CH4,molecular_ij [m^2/s],2.220814e-11

examples/methane_hydrogen_porous_media/run.py

@@ -0,0 +1,20 @@
+"""Calculate diffusion coefficients for methane (CH4)/H2 in porous media"""
+
+
+def main():
+    from porous_media.parameters import FluidMixture
+    void_fraction, d_pore, tortuosity = 0.35, 0.5e-9, 2.
+    names = ['CH4', 'H2']  # must match what in parameter file!
+    molecular_weights = [12.011+1.008*4, 1.008*2]
+    from raw_data.read_data import get_LJ_params, read_csv
+    LJ_data = read_csv('../../raw_data/LJparams.csv')
+    sigma, epsilon = get_LJ_params(names, LJ_data)
+    cls = FluidMixture(
+        names, molecular_weights, sigma, epsilon, void_fraction, d_pore, tortuosity
+    )
+    cls.write_params('input.csv')
+    cls.write_calculations('output.csv', 300., 101325.)
+
+
+if __name__ == '__main__':
+    main()

porous_media/__init__.py

@@ -0,0 +1,161 @@
+import math
+
+
+class PureFluid:
+    """Only one component is moving through porous media"""
+    def __init__(self, void_fraction, d_pore, tortuosity, molecular_weight=None):
+        """
+
+        :param molecular_weight:            molecular weight [kg/mol] for each component
+        :param void_fraction: void fraction of porous media
+        :param d_pore:        nominal pore diameter [meters]
+        :param tortuosity:    tortuosity of porous media
+        """
+        self.void_fraction = void_fraction
+        self.d_pore = d_pore
+        self.molecular_weight = molecular_weight
+        self.tortuosity = tortuosity
+        self.R = 8.314  # J/mol/K [=] m**3*Pa/mol/K
+
+    def knudsen(self, T):
+        """Knudsen diffusivity of component i; derived from kinetic theory of gases
+
+        Units:
+            R = m**3*Pa/mol/K [=] m*m*kg/mol/s/s/K
+
+            sqrt(RT/MW) [=] sqrt(m*m/s/s)
+
+
+        :param i: component name
+        :return:  Knudsen diffusivity m^2/s
+        """
+        return self.void_fraction*self.d_pore/self.tortuosity/3.*math.sqrt(8*self.R*T/math.pi/self.molecular_weight)
+
+    def write_params(self, file_name):
+        with open(file_name, 'w') as f:
+            for key, val in self.__dict__.items():
+                if isinstance(val, float):
+                    f.write('{},{},{}\n'.format(key, ' ', val))
+                elif isinstance(val, dict):
+                    f.write('{},{},{}\n'.format(key, ' ', ' '))
+                    for key2, val2 in val.items():
+                        if isinstance(key2, tuple):
+                            key2 = ' '.join(key2)
+                        f.write('{},{},{}\n'.format(' ', key2, val2))
+                elif isinstance(val, list):
+                    f.write('{},{},{}\n'.format(key, ' ', ' '))
+                    for val2 in val:
+                        f.write('{},{},{}\n'.format(' ', val2, ' '))
+                elif isinstance(val, str) or val is None:
+                    f.write('{},{},{}\n'.format(key, ' ', val))
+                else:
+                    raise Exception('Val type not found for {}:{}'.format(key, val))
+
+
+class FluidMixture(PureFluid):
+    """Multiple components flowing through porous media"""
+    def __init__(self, components, molecular_weight, sigma, epsilon_molecular, *args):
+        """
+
+        :param components:
+        :type components: list
+        :param molecular_weight:
+        :type molecular_weight: list
+        :param sigma:   collision diameter in Angstrom, get from Poling et al "The Properties of gases and liquids"
+        :type sigma: list
+        :param epsilon_molecular: molecular epsilon in K,get from Poling et al "The Properties of gases and liquids"
+        :type epsilon_molecular: list
+        :param args: args to pass to super
+        """
+        PureFluid.__init__(self, *args)
+        self.molecular_weight_i = {
+            key: val for key, val in zip(components, molecular_weight)
+        }
+        self.sigma_i = {
+            key: val for key, val in zip(components, sigma)
+        }
+        self.epsilon_i = {
+            key: val for key, val in zip(components, epsilon_molecular)
+        }
+        self.component_i = components
+
+    def knudsen_i(self, i, j, temperature, pressure):
+        """
+
+        :param i: component name
+        :return: Knudsen diffusivity for component *i*
+        """
+        self.molecular_weight = self.molecular_weight_i[i]
+        return self.knudsen(temperature)
+
+    def sigma_ij_rule(self, i, j):
+        """Lennard-Jones combining rule for sigma"""
+        return (self.sigma_i[i] + self.sigma_i[j]) / 2.
+
+    def epsilon_ij_rule(self, i, j):
+        """Lennard-Jones combining rule for epsilon"""
+        return math.sqrt(self.epsilon_i[i] * self.epsilon_i[j])
+
+    def omega_ij(self, w):
+        """Bird 1960 p. 746
+
+        :return: Temperature-dependent collision integral (dimensionless)
+        """
+
+        A1, B1, A2, B2 = 1.06036, 0.15610, 0.19300, 0.47635
+        A3, B3, A4, B4 = 1.03587, 1.52996, 1.76474, 3.89411
+        value = (
+                A1 / pow(w, B1)
+                + A2 / math.exp(B2 * w)
+                + A3 / math.exp(B3 * w)
+                + A4 / math.exp(B4 * w)
+        )
+        assert 0.5 < value < 2.7, 'Value predicted of %2.3f is not reasonable' % value
+        return value
+
+    def molecular_ij(self, i, j, T, P):
+        """Molecular diffusivites estimated by Chapman-Enskog
+
+        :param i: sorbate i name
+        :param j: sorbate j name
+        :param T: temperature in K
+        :param P: pressure in atm
+        :return: binary diffusion coefficient (m^2/s)
+        """
+
+        # convert molecular weights to g/mol to apply formula
+        M_i = self.molecular_weight_i[i]*1000.
+        M_j = self.molecular_weight_i[j]*1000.
+
+        return 1.858e-3 * math.sqrt(T*T*T * (1. / M_i + 1. / M_j)) / (
+                P * self.sigma_ij_rule(i, j) * self.sigma_ij_rule(i, j) *
+                self.omega_ij(T / self.epsilon_ij_rule(i, j))
+        )/100./100.
+
+    def effective_macropore_i(self, i, j, temperature, pressure):
+        """
+
+        :param temperature: temperature in K
+        :param pressure: pressure in Pa
+        :param i:       component i name
+        :param j:       other_component name
+        :return:        effective macropore diffusivity m^2/s
+        """
+        P_atm = pressure/101325.
+        return self.void_fraction/self.tortuosity/(
+                1. / self.molecular_ij(i, j, temperature, P_atm) + 1. / self.knudsen_i(i, j, temperature, pressure)
+        )
+
+    def write_calculations(self, output_file, temperature, pressure):
+        with open(output_file, 'w') as f:
+            for i in self.component_i:
+                for j in self.component_i:
+                    if i == j:
+                        continue
+                    for attr in [
+                        'effective_macropore_i',
+                        'knudsen_i',
+                        'molecular_ij'
+                    ]:
+                        func = getattr(self, attr)
+                        f.write('%s,%s,%s [m^2/s],%e\n' % (i, j, attr, func(i, j, temperature, pressure)))

porous_media/parameters.py

@@ -0,0 +1,16 @@
+"Substance","Name","sigma [angstrom]","epsilon [K]"
+Ar,"Argon",3.542,93.3
+He,"Helium",255.1,10.2
+Kr,"Krypton",3.655,178.9
+CH4,"Methane",3.758,148.6
+H2,"Hydrogen",2.827,59.7
+H2O,"Water",2.641,809.1
+N2,"Nitrogen",3.798,71.4
+CO2,"Carbon Dioxide",3.971,195.2
+O2,"Oxygen",3.467,106.7
+C2H6,"Ethane",4.443,215.7
+C2H4,"Ethylene",4.163,224.7
+C3H8,"Propane",5.118,237.1
+C3H6,"Propylene",4.678,298.9
+nC4H10,"Normal Butane",4.687,531.4
+iC4H10,"Iso Butane",5.278,330.1

raw_data/LJparams.csv

@@ -0,0 +1,3 @@
+LJparams
+========
+* found from Cussler "Diffusion", 3rd edition; parameters found from viscosities

raw_data/README.md

@@ -0,0 +1,27 @@
+def read_csv(f_name):
+    with open(f_name) as f:
+        header = list(map(eval, next(f).rstrip('\n').split(',')))
+        data = {key: [] for key in header}
+        for line in f:
+            for key, val in zip(header, line.rstrip('\n').split(',')):
+                try:
+                    data[key].append(float(val))
+                except ValueError:
+                    data[key].append(val)
+    return data
+
+
+def get_LJ_params(components, LJ_data):
+    """Get LJ params from csv file
+
+    :param components:
+    :return: LJ params sigma and epsilon
+    """
+    sigma = []
+    epsilon = []
+    for key in components:
+        assert key in LJ_data['Substance'], 'No parameters found for {} in LJparams.csv'.format(key)
+        index = LJ_data['Substance'].index(key)
+        sigma.append(LJ_data['sigma [angstrom]'][index])
+        epsilon.append(LJ_data['epsilon [K]'][index])
+    return sigma, epsilon