import numpy as np import scipy.linalg as lin import scipy.constants as const import pandas as pd import os, sys, platform sys.path.append(os.path.dirname(os.getcwd()) + '/soret_model/cpp/release') sys.path.append(os.path.dirname(os.getcwd()) + '/soret_model') sys.path.append(os.path.dirname(os.getcwd()) + '/soret_model/cpp') if (platform.system() == 'Linux'): from release_ubuntu.KineticGas import cpp_KineticGas else: from release_mac.KineticGas import cpp_KineticGas class KineticGas: computed_points = {} def __init__(self, comps, eos, mole_fracs=None, T=300, p=1e5, N=5): ''' :param comps (str): Comma-separated list of components, following Thermopack-convention :param eos (thermopack): An initialized equation of state, only used to get mole weights :param mole_fracs (array-like): list of mole fractions :param N (int > 0): Order of approximation. Be aware that N > 10 can be detrimental to runtime. This should be implemented in C++ or Fortran. ''' # Packing out variables in an n-component compatible way complist = comps.split(',') self.mole_weights = np.array([eos.compmoleweight(eos.getcompindex(comp)) for comp in complist]) self.mole_fracs = mole_fracs self.m0 = np.sum(self.mole_weights) self.M = self.mole_weights / self.m0 self.sigmaij = self.get_hard_sphere_radius(comps) self.sigma = np.diag(self.sigmaij) # Calculate the (temperature independent) soret-coefficient at the ideal gas state self.N = N binary_mole_frac_matr = np.zeros_like(self.sigmaij) for i in range(len(binary_mole_frac_matr)): for j in range(len(binary_mole_frac_matr)): if (mole_fracs[i] + mole_fracs[j] < 1e-10): binary_mole_frac_matr[i, j] = 0 else: binary_mole_frac_matr[i, j] = mole_fracs[i] / (mole_fracs[i] + mole_fracs[j]) binary_KineticGas_matr = np.empty_like(self.sigmaij, dtype=cpp_KineticGas) for i in range(1, len(binary_KineticGas_matr)): for j in range(i): binary_mole_weights = [self.mole_weights[i], self.mole_weights[j]] binary_sigmaij = [[self.sigmaij[i, i], self.sigmaij[i, j]], [self.sigmaij[j, i], self.sigmaij[j, j]]] binary_mole_fracs = [binary_mole_frac_matr[i, j], binary_mole_frac_matr[j, i]] binary_KineticGas_matr[i, j] = cpp_KineticGas(binary_mole_weights, binary_sigmaij, binary_mole_fracs, T, N) binary_alpha_T0_matr = np.zeros_like(self.sigmaij) for i in range(1, len(binary_KineticGas_matr)): for j in range(i): A = binary_KineticGas_matr[i, j].A_matrix A = np.array(A) delta = binary_KineticGas_matr[i, j].delta_vector try: d = lin.solve(A, delta) except (ValueError, np.linalg.LinAlgError): binary_alpha_T0_matr[i, j] = np.nan binary_alpha_T0_matr[j, i] = np.nan continue d_1, d0, d1 = d[N - 1], d[N], d[N + 1] binary_alpha_T0_matr[i, j] = - (5 / (2 * d0)) * ( (binary_mole_frac_matr[i, j] * d1 / np.sqrt(self.M[i])) + ( binary_mole_frac_matr[j, i] * d_1 / np.sqrt(self.M[j]))) binary_alpha_T0_matr[j, i] = - (binary_mole_frac_matr[i, j] / binary_mole_frac_matr[j, i]) * \ binary_alpha_T0_matr[i, j] self.binary_d_matr = np.zeros((len(self.sigmaij), len(self.sigmaij), 3)) for i in range(1, len(binary_KineticGas_matr)): for j in range(i): A = binary_KineticGas_matr[i, j].A_matrix A = np.array(A) delta = binary_KineticGas_matr[i, j].delta_vector try: d = lin.solve(A, delta) except (ValueError, np.linalg.LinAlgError): binary_alpha_T0_matr[i, j] = np.nan binary_alpha_T0_matr[j, i] = np.nan continue self.binary_d_matr[i, j, 0] = d[N - 1] self.binary_d_matr[i, j, 1] = d[N] self.binary_d_matr[i, j, 2] = d[N + 1] self.kT_vec = np.zeros_like(self.mole_weights) for i in range(len(binary_alpha_T0_matr)): for j in range(len(binary_alpha_T0_matr)): self.kT_vec[i] += (mole_fracs[i] + mole_fracs[j]) * binary_alpha_T0_matr[i, j] def alpha_T0(self, T): return self.kT_vec / (np.array(self.mole_fracs) * (1 - np.array(self.mole_fracs))) def get_hard_sphere_radius(self, comps): ''' Get hard-sphere diameters, assumed to be equal to Mie-potential sigma parameter. Gets Mie-parameters from the file mie_dev.xlsx. :param comps: Comma seperated list of components :return: N x N matrix of hard sphere diameters, where sigma_ij = 0.5 * (sigma_i + sigma_j), such that the diagonal is the radius of each component, and off-diagonals are the average diameter of component i and j. ''' df = pd.read_excel(os.path.dirname(__file__) + '/mie_dev.xlsx') sigma_i = np.array([df.loc[df['comp'] == comp]['sigma'].iloc[0] for comp in comps.split(',')]) epsilon = np.array([df.loc[df['comp'] == comp]['epsilon'].iloc[0] for comp in comps.split(',')]) epsilon *= const.Boltzmann sigma_i *= 1e-10 sigma = [quad(self.BH_integrand, 0, sigma_i[i], args=(sigma_i[i], epsilon[i]))[0] for i in range(len(sigma_i))] sigma_i = sigma sigma_ij = 0.5 * np.sum(np.meshgrid(sigma_i, np.vstack(sigma_i)), axis=0) return sigma_ij def get_sigma(self, comps): df = pd.read_excel(os.path.dirname(__file__) + '/mie_dev.xlsx') sigma_i = np.array([df.loc[df['comp'] == comp]['sigma'].iloc[0] for comp in comps.split(',')]) * 1e-10 sigma_ij = 0.5 * np.sum(np.meshgrid(sigma_i, np.vstack(sigma_i)), axis=0) return sigma_ij def BH_integrand(self, r, sigma, epsilon): lambda_r = 12 lambda_a = 6 return (1 - np.exp(-self.u_Mie(r, sigma, epsilon, lambda_r, lambda_a) / (self.T * const.Boltzmann))) def u_Mie(self, r, sigma, epsilon, lambda_r, lambda_a): C = lambda_r / (lambda_r - lambda_a) * (lambda_r / lambda_a) ** (lambda_a / (lambda_r - lambda_a)) return C * epsilon * ((sigma / r) ** lambda_r - (sigma / r) ** lambda_a)