Berechnung von Impuls- und Wärmeübertragungseigenschaften einfacher Gase nach kinetischer Gastheorie

lvkg.py

import numpy as np
from scipy.interpolate import interp1d
from scipy import interpolate
t = 25 + 273.15 # K
p = 1.01325 # bar
komponente = np.array(['N2','O2','Ar','Luft'])
# Bird Tabelle E.1
tc = np.array([126.2, 154.4, 150.7, 132]) # K
pc = np.array([33.5, 49.7, 48.0, 36.4]) # atm
muc = np.array([180., 250., 264., 193.])*1e-6 # g/cm/s
# Pa s= N/m^2 s=kg m/s^2/m^2 s = 
# kg/m/s * 1000g/g * 1m/100cm = 10 g/cm/s
kc = np.array([86.8,105.3,71.0,90.8])*1e-6*4.184*100 # cal/cm/s/K
# * 4.184J/cal * 100cm/m [=] W/m/K
tr = t/tc
pr = p/pc
# Bird Fig. 1.3-1
mur = np.array([0.97,0.85,0.86,0.94])
kr = np.array([0.675,0.58, 0.595,0.65])
mu = mur * muc
k = kr * kc
y_i = np.array([78, 21, 1, 0])/100.
mm_i = np.array([28, 32, 40, 28.97])
mm = sum(y_i * mm_i)
def phi_alpha_beta(mm_i, mu):
    phi_ab = np.zeros([mm_i.size, mu.size])
    for alpha in range(phi_ab.shape[0]):
        for beta in range(phi_ab.shape[1]):
            phi_ab[alpha, beta] = 1/np.sqrt(8)*(
                1+mm_i[alpha]/mm_i[beta])**(-1/2.)*(
                1+(mu[alpha]/mu[beta])**(1/2.)*(
                    mm_i[beta]/mm_i[alpha]
                )**(1/4.)
            )**2
    return phi_ab
rho = interpolate.interp1d(
    [10,20,30],
    [1.231, 1.189, 1.149]
)(25)
cp = interpolate.interp1d(
    [10,20,30],
    [1.006,1.006,1.007]
)(25)
cv = interpolate.interp1d(
    [10,20,30],
    [0.7174,0.7178,0.7183]
)(25)
lambda_g = interpolate.interp1d(
    [10,20,30],
    [25.12,25.87,26.62]
)(25)
eta = interpolate.interp1d(
    [10,20,30],
    [x*1e-6 for x in [17.72, 18.21, 18.69]]
)(25)
r_berechnet = (1.0065-0.71805)*28.96
print('Luft Zusammensetzung: ')
print('N2:O2:Ar:(Luft), ' + str(y_i*100))
print('Molare Masse, mm/(g/mol)')
print(mm)
print('Reduzierte Eigenschaften bei T, P')
print('Tr: ' + str(tr))
print('Pr: ' + str(pr))
print('mur: ' + str(mur) + '(Bird Fig. 1.3-1)')
print('mu, g/cm/s')
print(mu)
print('kr: ' + str(kr) + '(Bird Fig. 9.2-1)')
print('k, W/m/s')
print(k)
print('Stoffdaten (VDI Wärmeatlas) bei 1atm, 25°C')
print('Dichte, rho/(kg / m^3)')
print(rho)
print('cp kJ/(kg K)')
print(cp)
print('cv kJ/(kg K)')
print(cv)
print('R=cp-cv (für monatomige Gase)')
print(r_berechnet)
print('Wärmeleitfähigkeit, lambda/(mW/(m K))')
print(lambda_g)
print('Dynamische Viscosität, direkt abgelesen \n'+
      '(VDI-Wärmeatlas) eta/(Pa s)')
print('{:0.4e}'.format(eta.item()))
print('Phi_alpha_beta, nach Bird Gl. 1.4-16')
print(phi_alpha_beta(mm_i,mu))
print('Dynamische Viskosität, nach dem Prinzip der \n'+
      'korrespondierenden Zuständen, mu/(g/cm/s)')
print('{:0.4e}'.format(
    sum(
        y_i * mu / phi_alpha_beta(mm_i,mu).dot(y_i)
    ).item())
)
print('Wärmeleitfähigkeit, nach dem Prinzip der \n'+
      'korrespondierenden Zuständen, k/(W/m/K)')
print('{:0.4e}'.format(
    sum(
        y_i * k / phi_alpha_beta(mm_i,k).dot(y_i)
    ).item())
)
print('Dynamische Viskosität, bei gegebenen Stoffdaten \n'+
      'der reinen Stoffe (Viskositäten aus VDW Wärmeatlas) \n' +
      ', mu/(g/cm/s)')
mu = np.array([
            interpolate.interp1d(
                [20,30],
                [17.57e-6*10,18.03e-6*10]
            )(25),
            interpolate.interp1d(
                [20,30],
                [20.27e-6*10,20.55e-6*10]
            )(25),
            interpolate.interp1d(
                [0,25],
                [21.1e-6*10,31.6e-6*10]
            )(25),
            interpolate.interp1d(
                [0,25],
                [17.22e-6*10,18.45e-6*10]
            )(25)
        ])
print('{:0.4e}'.format(sum(
    y_i * mu / phi_alpha_beta(mm_i,mu
    ).dot(y_i))))
print('1 Pa s = 1000/100 g/cm/s')
print('Transporteigenschaften nach kinetischer Gastheorie:\n'+
      '1. Maxwell 1860, Stoßzylinder mit harten Kugeln und\n'+
      ' mittlerer freien Weglänge.')
c = np.sqrt(8*1.3806e-23*t/(np.pi*mm_i/1000./6.0221e23))
# sqrt(kg m^2/s^2/K*K/kg) = m/s
sigma = np.array([0.43, 0.40, 0.36,np.nan])*1e-9**2
lambda_mfw = 1.3806e-23*t/(np.sqrt(2)*sigma*1)*1/101325.
cv=np.array([0.743,0.6584,0.3119, np.nan]) # kJ/(kg K)
mu=1/3.*lambda_mfw*c*101325./(8.3145*t)*mm_i/1000.
k=1/3.*lambda_mfw*c*101325./(8.3145*t)*cv*mm_i
print('Mittlere Geschwindigkeit, nach Bird Gl. 1.4-1')
print('u, m/s')
print(c)
print('Mittlere freie Weglänge, nach Bird Gl. 1.4-3')
print('lambda, m')
print(lambda_mfw)
print('Stoßquerschnitt (Atkins T. 21.1)')
print('sigma, m^2')
print(sigma)
print('Wärmekapazität bei konstantem Volumen (VDI Wärmeatlas)')
print('Cv, J/g/K')
print(cv)
print('Dynamische Viskosität, mu/(g/cm/s)')
print('{:0.4e}'.format(sum(
    y_i[:-1] * mu[:-1] / phi_alpha_beta(mm_i[:-1],mu[:-1]
    ).dot(y_i[:-1]))))
print('Wärmeleitfähigkeit, mW/m/K')
print('{:0.4f}'.format(sum(
    y_i[:-1] * k[:-1]*1000 / phi_alpha_beta(mm_i[:-1],k[:-1]
    ).dot(y_i[:-1]))))
print('Ergebnisse weisen eine wesentliche Abweichung dar, \n'+
      'vermutlich wegen schwacher Temperaturabhängigkeit \n'+
      '(Wurzel T), und wegen Kugelförmiger Annahme ohne \n'+
      'Wechselwirkungen.'
     )
print('')
print('2. Chapman Enskog, monatomige Gasen unter \n'+
      'Berücksichtigung der Lennard-Jones-Parameter\n'+
      '(anziehende und abstoßende Wechselwirkungen)\n'+ 
      'und Stoßfunktion.')
# Lennard-Jones Parameter (Bird Tabelle E.1)
l_j_epsilon_d_k = np.array([99.8,113,122.4,97.0]) # K
l_j_sigma =  np.array([3.667,3.433,3.432,3.617]) # Angstrom
k_t_d_e = t/l_j_epsilon_d_k
# Stoßintegral (Bird Tabelle E.2)
stossintegral_k_mu = interp1d(
    [2.40,2.50,2.60,2.7,2.8,2.9,3.0,3.1],
    [1.107,1.0933,1.0807,1.0691,1.0583,1.0482,1.0388,1.0300]
)(k_t_d_e)
konst_1 = 5/16*np.sqrt(
    1.3806e-23*1000*100**2/6.022e23/np.pi
)*(10**10/100)**2 # 1/cm/s
mu = konst_1 * np.sqrt(mm_i*t)/(
    l_j_sigma**2*stossintegral_k_mu)
konst_2=(9/4*8.3145+cv*mm_i)*1/4.184*konst_1
k = konst_2 * np.sqrt(t/mm_i)/(
    l_j_sigma**2*stossintegral_k_mu)*4.184*100
print('Lennard-Jones Parameter (Bird Tabelle E.1)')
print('epsilon/k, 1/K')
print(l_j_epsilon_d_k)
print('sigma, Angstrom')
print(l_j_sigma)
print('kT/epsilon')
print(k_t_d_e)
print('Stoßintegral, Omega_k=Omega_mu')
print(stossintegral_k_mu)
print('Dynamische Viskosität, nach Bird Gl. 1.4-14')
print('mu/(g/cm/s)')
print('{:0.4e}'.format(sum(
    y_i * mu / phi_alpha_beta(mm_i,mu
    ).dot(y_i))))
print('Wärmeleitfähigkeit, nach Bird Gl. 9.3-13')
print('k,  mW/m/K')
print('{:0.4f}'.format(sum(
    y_i[:-1] * k[:-1]*1000 / phi_alpha_beta(mm_i[:-1],k[:-1]
    ).dot(y_i[:-1]))))

phialphabeta.py

print('Bird E 1.4-2')
komponente=np.array(['CO2','O2','N2'])
mm_i=np.array([44.01,32,28.02])
mu=np.array([1462e-7,2031e-7,1754e-7])
x_a=np.array([0.133,0.039,0.828])
phi_ab = np.zeros([komponente.size, komponente.size])
for alpha in range(phi_ab.shape[0]):
    for beta in range(phi_ab.shape[1]):
        phi_ab[alpha, beta] = 1/np.sqrt(8)*(
            1+mm_i[alpha]/mm_i[beta])**(-1/2.)*(
            1+(mu[alpha]/mu[beta])**(1/2.)*(
                mm_i[beta]/mm_i[alpha]
            )**(1/4.)
        )**2
print('phi_alpha_beta')
print(phi_ab)
print('sum(x_alpha phi_alpha_beta)')
print(phi_ab.dot(x_a))
print('Viskosität, g/cm/s')
print('{:g}'.format(
    sum(x_a * mu / phi_ab.dot(x_a))*1e7
)+'(10)^-7 g/cm/s')