# Properties of gas mixture at 1000K

## Viscosity 

Pure monatomic gas [1] with mass of a single molecule $m=M/N_A$ where $M$ is the molar mass and $N_A$ the Avogadro constant.

$\begin{align}
\mu=&\frac{5}{16}\frac{\sqrt{\pi m \kappa T}}{\pi\sigma^2\Omega_\mu}\\
=&\frac{5}{16}\sqrt{\frac{1.3806\cdot 10^{-23}}{\pi\cdot 1000\cdot 6.022\cdot 10^{23}}}\frac{1}{(10^{-10})^2}\frac{kg}{m\cdot s}\cdot\frac{\sqrt{\left(\frac{M}{g/mol}\right)\left(\frac{T}{K}\right)}}{\left(\frac{\sigma}{\unicode{x212B}}\right)^2\Omega_\mu}\\
=&\underset{2.6696\cdot 10^{-5}}{\underbrace{\frac{5}{16}\sqrt{\frac{1.3806\cdot 10^{-23}}{\pi\cdot 1000\cdot 6.022\cdot 10^{23}}}\frac{1}{(10^{-10})^2}\frac{1000}{100}}}\frac{g}{cm\cdot s}\cdot\frac{\sqrt{\left(\frac{M}{g/mol}\right)\left(\frac{T}{K}\right)}}{\left(\frac{\sigma}{\unicode{x212B}}\right)^2\Omega_\mu}
\end{align}$

with $\kappa$ Boltzmann constant, $\sigma$ Lennard-Jones collision diameter, $\epsilon$ Lennard-Jones maximum energy of atraction between a pair of molecules, $\Omega_\mu$ collision integral for viscosity.

Empirical equation for the collision integral for viscosity [2]:

$\Omega_\mu=A\cdot(T^*)^{-B}+C\cdot exp(-D\cdot(T^*))+E\cdot exp(-F\cdot(T^*))$

$(T^*)=\kappa T/\epsilon, A=1.16145, B=0.14874, C=0.52487, D=0.77320, E=2.16178, and F=2.43787$ applicable for $0.3\le(T^*)\le 100$.

Viscosity of the gas mixture with molar fractions $x_i$ at low density by the semiempirical relation [1]:

$\mu_{mix}=\sum_{\alpha=1}^{N}\frac{x_\alpha \mu_\alpha}{\sum_\beta x_\beta\Phi_{\alpha\beta}}$

$\Phi_{\alpha\beta}=\frac{1}{\sqrt{8}}\left(1+\frac{M_\alpha}{M_\beta}\right)^{-1/2}\left[1+\left(\frac{\mu_\alpha}{\mu_\beta}\right)^{1/2}\left(\frac{M_\beta}{M_\alpha}\right)^{1/4}\right]^2$

Lennard-Jones parameters [2]

i | $\sigma/\unicode{x212B}$ | $(\frac{\epsilon}{\kappa}) / K$
- | - | -
$H_2$ | 2.827 | 59.7
$CH_4$ | 3.758 | 148.6
$CO_2$ | 3.941 | 195.2
$CO$ | 3.690 | 91.7
$H_2O$ | 2.641 | 809.1
$N_2$ | 3.798 | 71.4
$O_2$ | 3.467 | 106.7

## Density

Ideal gas.

$\rho=\frac{M}{V}=\frac{p}{RT}M$

Density of mixture with molar fractions $x_i$:

$\rho=\frac{p}{RT}\sum_i x_i M_i$

## Heat capacity (constant pressure)

"NASA-polynomials" with coefficients from [3], applicable between 200K and 1000K, and between 1000K and 6000K.

$\frac{Cp}{R}=a_1+a_2 T+a_3 T^2+a_4 T^3+a_5 T^4$

Heat capacity of the gas mixture with molar fractions $x_i$:

$Cp_{mix}=\sum_i x_i Cp_i$

## Thermal conductivity

Three correlations to test, based on Eucken's semiempirical relation for polyatomic gases at low density, see ref. [1] and [2]. Using eq. 10-3.2 in [2] as basis:

$\frac{\lambda M}{\eta C_v}=\frac{15/4}{C_v/R}+f_{int}\left(1-\frac{3/2}{C_v/R}\right)$

Here, $M$ is the molar mass, $\mu$ is the dynamic viscosity, $C_v$ is the constant volume heat capacity, $\frac{\lambda M}{\eta C_v}$ is the Eucken number and $f_{int}$ represents the contribution due to internal degrees of freedom. The contribution due to translational degrees of freedom, $f_{tr}$, is set to $f_{tr}=2.5$, because the limiting form for monatomic gases at low pressure is $\frac{\lambda M}{\eta C_v}=2.5$. Additionally, the relationship between $C_v$ and $C_p$ for an ideal gas is: $C_p=C_v-T \left(\frac{\partial V}{\partial T}\right)^2_P/\left(\frac{\partial V}{\partial P}\right)_T=C_v+R$.

0. Monatomic gas: $C_v=\frac{3}{2}R$ and $f_{int}=0$

$\frac{\lambda M}{\eta C_v}=\frac{15/4}{3/2R}R=\frac{5}{2}$

1. Polyatomic gas by Eucken correlation: $f_{int}=1$

$\begin{align}
\lambda&=\left(\frac{15}{4}R+C_v-\frac{3}{2}R+\right)\frac{\mu}{M}=\left(\frac{9}{4}R+C_p-R\right)\frac{\mu}{M}\\
&=\left(Cp+\frac{5}{4}R\right)\mu\\
&=\left(\frac{Cp}{J/mol/K}+\frac{5}{4}\frac{R}{J/mol/K}\right)\frac{\mu}{kg/m/s}\frac{g/mol}{M}\cdot \frac{J}{mol K}\frac{kg}{m s}
\frac{mol}{g}\frac{1000g}{1kg}[=]\frac{W}{m K}\end{align}$

2. Polyatomic gas by modified Eucken correlation: $f_{int}\approx 1.32$

$\begin{align}
\lambda&=\left(1.32 C_p+\left[\frac{15}{4}-1.32-1.32\cdot\frac{3}{2}\right]R\right)\frac{\mu}{M}\\
&=\left(1.32\frac{Cp}{J/mol/K}+0.45\frac{R}{J/mol/K}\right)\frac{\mu}{kg/m/s}\frac{g/mol}{M}\cdot \frac{J}{mol K}\frac{kg}{m s}
\frac{mol}{g}\frac{1000g}{1kg}[=]\frac{W}{m K}\end{align}$

3. Polyatomic non-polar gas compromise (Stiel, Thodos [2]): $f_{int}\approx 1.15$

$\begin{align}
\lambda&=\left(1.15 C_p+\left[\frac{15}{4}-1.15-1.15\cdot\frac{3}{2}\right]R\right)\frac{\mu}{M}\\
&=\left(1.15\frac{Cp}{J/mol/K}+0.875\frac{R}{J/mol/K}\right)\frac{\mu}{kg/m/s}\frac{g/mol}{M}\cdot \frac{J}{mol K}\frac{kg}{m s}
\frac{mol}{g}\frac{1000g}{1kg}[=]\frac{W}{m K}\end{align}$

Conductivity of the gas mixture with molar fractions $x_i$ at low density by the semiempirical relation[1]:

$\lambda_{mix}=\sum_{\alpha=1}^{N}\frac{x_\alpha \lambda_\alpha}{\sum_\beta x_\beta\Phi_{\alpha\beta}}$

$\Phi_{\alpha\beta}=\frac{1}{\sqrt{8}}\left(1+\frac{M_\alpha}{M_\beta}\right)^{-1/2}\left[1+\left(\frac{\mu_\alpha}{\mu_\beta}\right)^{1/2}\left(\frac{M_\beta}{M_\alpha}\right)^{1/4}\right]^2$

Example with two components

$\lambda_{mix}=\frac{x_1\lambda_1}{x_1+x_2\frac{1}{\sqrt{8}}\left(1+\frac{M_1}{M_2}\right)^{-1/2}\left[1+\left(\frac{\mu_1}{\mu_2}\right)^{1/2}\left(\frac{M_2}{M_1}\right)^{1/4}\right]^2}+\frac{x_2\lambda_2}{x_1\frac{1}{\sqrt{8}}\left(1+\frac{M_2}{M_1}\right)^{-1/2}\left[1+\left(\frac{\mu_2}{\mu_1}\right)^{1/2}\left(\frac{M_1}{M_2}\right)^{1/4}\right]^2+
x_2}$

since $\frac{1}{\sqrt{8}}\left(1+1\right)^{-1/2}\left[1+1\right]^2=1$

## Diffusion coefficient

Binary system, gas A diffusing through gas B at low to moderate pressure, with mean free path $\lambda=\frac{kT}{\sqrt{2}(\pi\sigma_{AB})}$, collision cross section $(\pi\sigma_{AB})$, mean velocity $\bar{c}=\sqrt{\frac{8RT}{\pi \tilde{M}_{AB}}}$ and number density of molecules in the mixture $\mathcal{N}=\frac{n}{V}N_A=\frac{p}{RT}N_A$. The characteristic length of the binary system for these purposes is given by the mean of both collision cross section diameters $\sigma_{AB}=\frac{1}{2}(\sigma_A+\sigma_B)$, the depth of the potential energy minimum by the geometric mean $\epsilon=(\epsilon_A\cdot\epsilon_B)^{1/2}$ and the molecular mass by $\frac{1}{M_{AB}}=\frac{1}{2}\cdot(\frac{1}{M_{A}}+\frac{1}{M_{B}})$. Replacing the Boltztmann constant $k=R/N_A$ into the number density of molecules, $\mathcal{N}=\frac{p}{k T}$; the molar mass by the molecular mass $M_{AB}=\tilde{M}_{AB}/N_A$, and a first approach to $D_{AB}$ from the kinetic gas theory yields eq. (1), which has comparable form to eq. 11.3-1 in ref. [2].

$\begin{align}
D_{AB}&=\frac{1}{3}\lambda \bar{c}=\frac{1}{3}\frac{kT}{\sqrt{2}(\pi\sigma_{AB})}\cdot\sqrt{\frac{8RT}{\pi \tilde M_{AB}}}=\frac{1}{3\pi}\frac{\sqrt{4\pi k T/M_{AB}}}{\mathcal{N}\cdot\pi\sigma_{AB}^2}\tag{1}
\end{align}$

Substituting constant terms to display functionality with pressure, temperature and molar mass,

$\begin{align}
D_{AB}&=\underbrace{\frac{1}{3\pi}\sqrt{\frac{4}{\pi}}}_{0.1197}\cdot k^{3/2}N_A^{1/2}\cdot\frac{T^{3/2}}{p\cdot \tilde{M}_{AB}^{1/2}\cdot \sigma_{AB}^2}\\
&=\frac{1}{3\pi}\sqrt{\frac{4}{\pi}}\cdot k^{3/2}N_A^{1/2}\cdot
\frac{\left(\frac{T}{K}\right)^{3/2}}{\left(\frac{p}{bar}\right)\cdot \left(\frac{\tilde{M}_{AB}}{g/mol}\right)^{1/2}\cdot \left(\frac{\sigma_{AB}}{\mathring A}\right)^2}
\cdot\frac{K^{3/2}\cdot \left(\frac{1}{10^{-3}kg/mol}\right)^{1/2}}{10^5 Pa\cdot \left(10^{-10}m\right)^2}\cdot\left(\frac{100 cm}{m}\right)^2\\
&=\frac{1}{3\pi}\cdot\sqrt{\frac{4}{\pi}}\cdot\left(1.3806\cdot 10^{-23}\frac{J}{K}\right)^{3/2}\left(6.0221\cdot 10^{23}\frac{1}{mol}\right)^{1/2}\cdot
\frac{K^{3/2}\cdot \left(\frac{1}{10^{-3}kg/mol}\right)^{1/2}}{10^5 Pa\cdot \left(10^{-10}m\right)^2}\cdot\left(\frac{100 cm}{m}\right)^2\cdot
\frac{\left(\frac{T}{K}\right)^{3/2}}{\left(\frac{p}{bar}\right)\cdot \left(\frac{\tilde{M}_{AB}}{g/mol}\right)^{1/2}\cdot \left(\frac{\sigma_{AB}}{\mathring A}\right)^2}\\
&=\frac{\frac{1}{3\pi}\cdot\sqrt{\frac{4}{\pi}}\cdot\left(1.3806\cdot 10^{-23}\right)^{3/2}\left(6.0221\cdot 10^{23}\right)^{1/2}\left(10^3\right)^{1/2}\cdot 100^2}{10^5\cdot 10^{-20}} \cdot
\frac{\frac{kg^{3/2}m^3}{s^3}\cdot mol^{-1/2}\frac{mol^{1/2}}{kg^{1/2}}\frac{cm^{2}}{m^{2}}}{\frac{kg~m}{m^2 s^2}m^2}\cdot
\frac{\left(\frac{T}{K}\right)^{3/2}}{\left(\frac{p}{bar}\right)\cdot \left(\frac{\tilde{M}_{AB}}{g/mol}\right)^{1/2}\cdot \left(\frac{\sigma_{AB}}{\mathring A}\right)^2}\\
&=0.00150716\cdot\frac{\left(\frac{T}{K}\right)^{3/2}}{\left(\frac{p}{bar}\right)\cdot \left(\frac{\tilde{M}_{AB}}{g/mol}\right)^{1/2}\cdot \left(\frac{\sigma_{AB}}{\mathring A}\right)^2}~\frac{cm^2}{s}
\end{align}$

The Chapman-Enskog version, eq. (2) in ref. [2] has a similar form except with the factor $\frac{3}{16}$ instead of $\frac{1}{3\pi}$ (77 % higher), and adds the collision integral $\Omega_D$ as well as a correction term $f_D$, which lies between 1.0 and 1.02. For present purposes, $f_D=1.0$.

$\begin{align}
D_{AB}&=\frac{3}{16}\frac{\sqrt{4\pi k T/M_{AB}}}{\mathcal{N}\cdot\pi\sigma_{AB}^2\cdot\Omega_D}\cdot f_D\tag{2}\\
&=0.00266337\cdot\frac{\left(\frac{T}{K}\right)^{3/2}}{\left(\frac{p}{bar}\right)\cdot \left(\frac{\tilde{M}_{AB}}{g/mol}\right)^{1/2}\cdot \left(\frac{\sigma_{AB}}{\mathring A}\right)^2\cdot\Omega_D}~\frac{cm^2}{s}
\end{align}$

Empirical equation for the collision integral for diffusion [2]:

$\begin{align}
\Omega_D&=A\cdot(T^*)^{-B}+C\cdot exp(-D\cdot(T^*))+E\cdot exp(-F\cdot(T^*))+G\cdot exp(-H\cdot(T^*))\\
(T^*)&=\kappa T/\epsilon
\end{align}$

$\begin{array}{cc}
A=1.06036 & B=0.15610\\
C=0.19300 & D=0.47635\\
E=1.03587 & F=1.52996\\
G=1.76474 & H=3.89411\\
\end{array}$

For (dilute) component $i$ diffusing through a mixture of gases, Blanc's law [2] yields eq. (3) for the diffusion coefficient of component $i$ in the mixture, $D_{im}$, where $D_{ij}=D_{ji}$.

$\begin{align}
D_{im}=\frac{1}{\sum_{j=1\\i\neq j}^{n}{\frac{x_j}{D_{ij}}}}\tag{3}
\end{align}$

Eq. (4) is an equivalent form for somewhat less diluted gases ($x_i$ not quite 0) as found in [19] eq. 12.2.

$\begin{align}
D_{im}=\frac{1-x_i}{\sum_{j=1\\i\neq j}^{n}{\frac{x_j}{D_{ij}}}}\tag{4}
\end{align}$


```python
from numpy import pi,array,float64,exp,sqrt,piecewise,outer,linspace,diag,arange,ones,zeros,finfo
from matplotlib import cycler,pyplot as plt

R=8.31446261815324 # J/mol/K
k_B=1.380649e-23 # J/K
N_A=6.02214076e+23 # 1/mol
R=k_B*N_A

labels=['H2','CH4','CO2','CO','H2O','N2','O2','Ar','He','C6H6','C6H12','C6H14'] # order of components
sigma_lj=array([2.827,3.758,3.941,3.69,2.641,3.798,3.467,3.542,2.551,5.349,6.182,5.949]) # Lennard Jones sigma, Angström 
epsilon_lj_ov_k=array([59.7,148.6,195.2,91.7,809.1,71.4,106.7,93.3,10.22,412.3,297.1,399.3]) # Lennard Jones epsilon/Boltztmann k, K
M=array([2,16,44,28,18,28,32,40,4,78,84,86]) # g/mol
# coefficients of NASA polynomials from Burcat, A., & Ruscic, B. (2001). Third millennium ideal gas and condensed phase thermochemical database for combustion. Technion-Israel Institute of Technology.
a1_a7_low=array([x.split('\t') for x in """2,34433112	0,007980521	-1,95E-05	2,02E-08	-7,38E-12	-917,935173	0,683010238
5,14825732	-0,013700241	4,94E-05	-4,92E-08	1,70E-11	-10245,3222	-4,63322726
2,356813	0,00898413	-7,12E-06	2,46E-09	-1,43E-13	-48371,971	9,9009035
3,5795335	-0,000610354	1,02E-06	9,07E-10	-9,04E-13	-14344,086	3,5084093
4,1986352	-0,002036402	6,52E-06	-5,49E-09	1,77E-12	-30293,726	-0,84900901
3,53100528	-0,000123661	-5,03E-07	2,44E-09	-1,41E-12	-1046,97628	2,96747038
3,78245636	-0,002996734	9,85E-06	-9,68E-09	3,24E-12	-1063,94356	3,65767573
2,5	0	0,00E+00	0,00E+00	0,00E+00	-745,375	4,37967491
2,5	0	0,00E+00	0,00E+00	0,00E+00	-745,375	0,928723974
0,504818632	0,018502064	7,38E-05	-1,18E-07	5,07E-11	8552,47913	21,6412893
4,04357527	-0,006196083	1,77E-04	-2,23E-07	8,64E-11	-16920,3544	8,52527441
9,87121167	-0,00936699	1,70E-04	-2,15E-07	8,45E-11	-23718,5495	-12,4999353
""".replace(',','.').split('\n') if len(x)>0],dtype=float64) 
a1_a7_high=array([x.split('\t') for x in """2,93286575	0,000826608	-1,46E-07	1,54E-11	-6,89E-16	-813,065581	-1,02432865
1,911786	0,00960268	-3,38E-06	5,39E-10	-3,19E-14	-10099,2136	8,48241861
4,6365111	0,002741457	-9,96E-07	1,60E-10	-9,16E-15	-49024,904	-1,9348955
3,0484859	0,001351728	-4,86E-07	7,89E-11	-4,70E-15	-14266,117	6,0170977
2,6770389	0,002973182	-7,74E-07	9,44E-11	-4,27E-15	-29885,894	6,88255
2,95257637	0,0013969	-4,93E-07	7,86E-11	-4,61E-15	-923,948688	5,87188762
3,66096065	0,000656366	-1,41E-07	2,06E-11	-1,30E-15	-1215,97718	3,41536279
2,5	0	0,00E+00	0,00E+00	0,00E+00	-745,375	4,37967491
2,5	0	0,00E+00	0,00E+00	0,00E+00	-745,375	0,928723974
11,0809576	0,020717675	-7,52E-06	1,22E-09	-7,36E-14	4306,41035	-40,041331
13,214597	0,035824343	-1,32E-05	2,17E-09	-1,32E-13	-22809,2102	-55,3518322
19,5158086	0,026775394	-7,50E-06	1,20E-09	-7,52E-14	-29436,2466	-77,4895497
""".replace(',','.').split('\n') if len(x)>0],dtype=float64)

Cp_R_coefs_200_1000_K=a1_a7_low[:,:4+1] # for Cp function, coefficients from a1 to a5 are applicable
Cp_R_coefs_1000_6000_K=a1_a7_high[:,:4+1] # for Cp function, coefficients from a1 to a5 are applicable

def mu(T):
    T_=outer(T,1/epsilon_lj_ov_k).squeeze() # outer handles T as array. squeeze eliminates additional dimension when T has length 1
    MT=outer(T,M).squeeze() # outer handles T as array. squeeze eliminates additional dimension when T has length 1
    omega_mu=1.16145*T_**(-0.14874)+0.52487*exp(-0.77320*T_)+2.16178*exp(-2.43787*T_)
    return 5/16*sqrt(k_B/(pi*1000*N_A))*1/1e-10**2*sqrt(MT)/(sigma_lj**2*omega_mu)  # Pa s

def D_ij(T,p):
    """Returns binary diffusion coefficients for p in Pa, T in K."""
    eps_lj_ov_k_ij=sqrt(outer(epsilon_lj_ov_k,epsilon_lj_ov_k))
    sigma_lj_ij=1/2*(sigma_lj.reshape([len(sigma_lj),1])+sigma_lj.reshape([1,len(sigma_lj)]))
    M_ij=1/(1/2*(1/M.reshape([1,len(M)])+1/M.reshape([len(M),1])))
    Tstar=T/eps_lj_ov_k_ij
    omega_D=1.06036*Tstar**(-0.15610)+0.19300*exp(-0.47635*Tstar)+1.03587*exp(-1.52996*Tstar)+1.76474*exp(-3.89411*Tstar)
    c1=3/(16)*sqrt(4/pi)*(1.3806e-23)**(3/2)*(6.0221e23)**(1/2)*1000**(1/2)*100**2/(1e5*1e-10**2) # 0.0026633 cm^2/s
    return c1*T**(3/2)/((p/1e5)*M_ij**(1/2)*sigma_lj_ij**2*omega_D)*100**(-2)  # m^2/s
    
def cp(T):
    result=((200<=T)&(T<=1000))*R*Cp_R_coefs_200_1000_K.dot(pow(T,array([[0],[1],[2],[3],[4]],dtype=float))).squeeze()+\
            ((1000<T)&(T<=6000))*R*Cp_R_coefs_1000_6000_K.dot(pow(T,array([[0],[1],[2],[3],[4]],dtype=float))).squeeze()
    return result.T # ensure row dimension is T

T=1000. # K
viscosity=mu(T)
density=101325/8.3145/T*M/1000
heat_capacity=cp(T).squeeze()
conductivity=(heat_capacity+5/4*R)*viscosity/M*1000

# mixture properties according to Bird R. B., Stweart W. E., Lightfoot E. N. (2002). Transport phenomena. 2nd ed. John Wiley & Sons. New York. S. 26, 276, 864
x=array([0,0,0,0,0,0.79,0.21,0,0,0,0,0])
phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
viscosity_air=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
conductivity_air=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
heat_capacity_air=heat_capacity.dot(x)
density_air=101325/8.3145/T*M.dot(x)/1000
mol_mass_air=M.dot(x)

x=array([0.05,0,0,0,0.95,0,0,0,0,0,0,0])
phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
viscosity_1=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
conductivity_1=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
heat_capacity_1=heat_capacity.dot(x)
density_1=101325/8.3145/T*M.dot(x)/1000
mol_mass_1=M.dot(x)

x_mass=array([0.2,0.5,0.05,0.05,0.2,0,0,0,0,0,0,0]) # ref. CH4 with mass fractions {0.2, 0.5, 0.05, 0.05, 0.2, 0.0} for {"Hydrogen", "Methane", "Carbondioxide", "Carbonmonoxide", "Water", "Nitrogen"}
x=x_mass/M/sum(x_mass/M)
phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
viscosity_2=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
conductivity_2=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
heat_capacity_2=heat_capacity.dot(x)
density_2=101325/8.3145/T*M.dot(x)/1000
mol_mass_2=M.dot(x)

print(f'viscosity, density, heat capacity, conductivity at {T}K')
print('\t\tviscosity/(Pa s)\tdensity/(kg/m^3)\tCp heat cap./(J/mol/K)\tconductivity/(W/m/K)\tmol mass/(kg/mol)')
for i in range(len(M)):
    print('{:s}\t\t{:0.16g}\t{:0.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}'.format(labels[i],viscosity[i],density[i],heat_capacity[i],conductivity[i],M[i]))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:16.16g}'.format('N2:O2 79:21',viscosity_air,density_air,heat_capacity_air,conductivity_air,mol_mass_air))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:16.16g}'.format('H2O:H2 95:5',viscosity_1,density_1,heat_capacity_1,conductivity_1,mol_mass_1))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:16.16g}'.format('H2:..:H2O '+':'.join(['{:d}'.format(int(round(x*100,0))) for x in x if round(x*100,0)!=0]),viscosity_2,density_2,heat_capacity_2,conductivity_2,mol_mass_2))

x_2=x # store value

print('\n')
print('Example 23 from ref. [4] Stephan, P. u. a. (2013). VDI-Wärmeatlas. 11. bearb. und erw. Aufl. Berlin, Heidelberg: Springer-Verlag. S. 652.')
print('Benzene:Argon 25:75 mol, 100.6 °C, 1 bar')
T=273.15+100.6 # K
p=1e5 # Pa
viscosity=mu(T)
density=1e5/8.3145/T*M/1000
heat_capacity=cp(T).squeeze()
conductivity=(heat_capacity+5/4*R)*viscosity/M*1000

x=array([0,0,0,0,0,0,0,0.75,0,0.25,0,0]) # ref. CH4 with mass fractions {0.2, 0.5, 0.05, 0.05, 0.2, 0.0} for {"Hydrogen", "Methane", "Carbondioxide", "Carbonmonoxide", "Water", "Nitrogen"}
phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
viscosity_3=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
conductivity_3=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
heat_capacity_3=heat_capacity.dot(x)
density_3=1e5/8.3145/T*M.dot(x)/1000
conductivity_3_2=(heat_capacity_3+5/4*R)*viscosity_3/sum(x*M)*1000
print('\t\tviscosity/(Pa s)\tdensity/(kg/m^3)\tCp heat cap./(J/mol/K)\tconductivity/(W/m/K)')
for i in range(len(M)):
    print('{:s}\t\t{:0.16g}\t{:0.16g}\t{:16.16g}\t{:16.16g}'.format(labels[i],viscosity[i],density[i],heat_capacity[i],conductivity[i]))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}'.format('N2:O2 79:21',viscosity_air,density_air,heat_capacity_air,conductivity_air))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}'.format('H2O:H2 95:5',viscosity_1,density_1,heat_capacity_1,conductivity_1))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}'.format(':'.join([labels[i] for i in range(len(labels)) if x_2[i]!=0])+' '+':'.join(['{:d}'.format(int(round(x_2[i]*100,0))) for i in range(len(labels)) if round(x_2[i]*100,0)!=0]),viscosity_2,density_2,heat_capacity_2,conductivity_2))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g}'.format(':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0]),viscosity_3,density_3,heat_capacity_3,conductivity_3))
print('{:s}\t{:0.16g}\t{:0.16g}\t{:0.16g}\t{:0.16g} {:s}'.format(':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0]),viscosity_3,density_3,heat_capacity_3,conductivity_3_2,'<<-- as (Cp+(5/4)R)mu/M'))
```

    viscosity, density, heat capacity, conductivity at 1000.0K
    		viscosity/(Pa s)	density/(kg/m^3)	Cp heat cap./(J/mol/K)	conductivity/(W/m/K)	mol mass/(kg/mol)
    H2		1.955982693264806e-05	0.02437308316795958	30.30498670043659	0.3980235536840263	               2
    CH4		2.724647778117298e-05	0.1949846653436767	71.90360841381916	0.1401434278292605	              16
    CO2		3.915118921174217e-05	0.5362078296951108	54.35948267311666	0.05761681285305009	              44
    CO		4.029426499361585e-05	0.3412231643514342	33.19282721770769	0.06272364377917297	              28
    H2O		3.582127599676216e-05	0.2193577485116363	41.25830244395893	0.102789909125832	              18
    N2		3.948253432565145e-05	0.3412231643514342	32.71195844300502	0.06078200327628414	              28
    O2		4.76940081694413e-05	0.3899693306873533	34.8851767183697	0.0674844207325538	              32
    Ar		5.213425236731442e-05	0.4874616633591917	20.7861565453831	0.04063765241563151	              40
    He		4.416966540192154e-05	0.04874616633591917	20.7861565453831	0.3442940923505748	               4
    C6H6		2.359575797510093e-05	0.9505502435504238	212.0760221493862	0.06729906474997142	              78
    C6H12		2.006278674556335e-05	1.023669493054303	 318.00734471702	0.07843604349518311	              84
    C6H14		2.022522824247268e-05	1.048042576222262	332.6136056435308	0.08066730780115443	              86
    N2:O2 79:21	4.117207251436728e-05	0.3514598592819772	33.1683342808316	0.06216273272687121	           28.84
    H2O:H2 95:5	3.57378708570406e-05	0.2096085152444524	40.7106366567828	0.1125694746118567	            17.2
    H2:..:H2O 69:22:1:1:8	2.79084126300121e-05	0.08388129187607421	40.31407281055314	0.2852935734654439	6.883108821155876
    
    
    Example 23 from ref. [4] Stephan, P. u. a. (2013). VDI-Wärmeatlas. 11. bearb. und erw. Aufl. Berlin, Heidelberg: Springer-Verlag. S. 652.
    Benzene:Argon 25:75 mol, 100.6 °C, 1 bar
    		viscosity/(Pa s)	density/(kg/m^3)	Cp heat cap./(J/mol/K)	conductivity/(W/m/K)
    H2		1.028132637698654e-05	0.06435949929114045	29.21473875768902	0.2036104469846543
    CH4		1.338210887722712e-05	0.5148759943291236	39.00714497024956	0.04131744787477274
    CO2		1.84570267008103e-05	1.41590898440509	40.28929767925152	0.02126013559557998
    CO		2.078472903641881e-05	0.9010329900759664	29.2969246507713	0.02946235557924449
    H2O		1.32397952810201e-05	0.5792354936202641	34.05782208069962	0.03269560115199573
    N2		2.064242982712038e-05	0.9010329900759664	29.2203428120627	0.02920418803407894
    O2		2.430579207725035e-05	1.029751988658247	29.9005159030836	0.03060524131504677
    Ar		2.685909715655945e-05	1.287189985822809	20.7861565453831	0.02093615243114621
    He		2.332600894950746e-05	0.1287189985822809	20.7861565453831	0.181821777601301
    C6H6		9.504577526556898e-06	2.510020472354478	104.4104591962687	0.01398921951534246
    C6H12		8.70649094210535e-06	2.703098970227899	137.1563725790117	0.01529330901610682
    C6H14		8.205274194357598e-06	2.767458469519039	170.7907185743312	0.01728677596168146
    N2:O2 79:21	4.117207251436728e-05	0.3514598592819772	33.1683342808316	0.06216273272687121
    H2O:H2 95:5	3.57378708570406e-05	0.2096085152444524	40.7106366567828	0.1125694746118567
    H2:CH4:CO2:CO:H2O 69:22:1:1:8	2.79084126300121e-05	0.08388129187607421	40.31407281055314	0.2852935734654439
    Ar:C6H6 75:25	1.844963387364072e-05	1.592897607455726	41.69223220810449	0.01712788854069524
    Ar:C6H6 75:25	1.844963387364072e-05	1.592897607455726	41.69223220810449	0.01941323047607249 <<-- as (Cp+(5/4)R)mu/M


## Error analysis
Following errors are presented compared to tables of experimental values at 100°C and 500°C in ref. [4] Stephan, P. u. a. (2013). VDI-Wärmeatlas. 11. bearb. und erw. Aufl. Berlin, Heidelberg: Springer-Verlag. S. 652.


```python
# tabulated values from [4] Stephan, P. u. a. (2013). VDI-Wärmeatlas. 11. bearb. und erw. Aufl. Berlin, Heidelberg: Springer-Verlag. S. 652.
# labels=['H2','CH4','CO2','CO','H2O','N2','O2','Ar','He','C6H6','C6H12','C6H14'] # order of components
visc_tab_100C=   [10.4,13.5,18.4,20.9,12.4,21.1,24.7,27.0,23.2,9.4,8.1,8.1] # uPa s
visc_tab_500C=   [17.1,23.4,33.3,34.5,28.6,35.0,41.4,46.5,38.6,19.3,15.4,16.1] # uPa s
visc_tab_100C=[x*1e-6 for x in visc_tab_100C] # Pa s
visc_tab_500C=[x*1e-6 for x in visc_tab_500C] # Pa s
conduct_tab_100C=[0.2149,0.0453,0.0225,0.0302,0.0239,0.0309,0.0324,0.0209,0.1793,0.0169,0.0203,0.0203]
conduct_tab_500C=[0.3678,0.1184,0.0536,0.0534,0.0669,0.0551,0.0597,0.0362,0.2983,0.0618,0.0735,0.0753]
visc_tab_100C_air=21.90e-6 # Pa s
visc_tab_500C_air=36.53e-6 # Pa s
conduct_tab_100C_air=31.62e-3 # W/m/K
conduct_tab_500C_air=55.80e-3 # W/m/K


def conductivity_fun(heat_capacity,viscosity,M,method):
    """Calculate conductivity by Eucken correlation, modified Eucken correlation and Stiel Thodos modification (nonpolar compromise)."""
    if method.lower()=='eucken':
        return (heat_capacity+5/4*R)*viscosity/M*1000
    elif method.lower()=='eucken_mod':
        return (1.32*heat_capacity+(1.77-1.32)*R)*viscosity/M*1000
    elif method.lower()=='nonpolar':
        return (1.15*heat_capacity+(2.025-1.15)*R)*viscosity/M*1000


for method in ['Eucken','Eucken_mod','nonpolar']: # loop through all three methods to check correlation

    T=273.15+100 # K
    p=1e5 # Pa
    viscosity=mu(T)
    density=1e5/8.3145/T*M/1000
    heat_capacity=cp(T).squeeze()
    conductivity=conductivity_fun(heat_capacity,viscosity,M,method)

    viscosity_100=viscosity;conductivity_100=conductivity # store for comparison

    # air properties 100°C
    x=array([0,0,0,0,0,0.7812,0.2096,0.0092,0,0,0,0]) # N2,O2,Ar
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_100C_air=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_100C_air=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_100C_air=heat_capacity.dot(x)
    density_100C_air=1e5/8.3145/T*M.dot(x)/1000
    conductivity_100C_air_2=conductivity_fun(heat_capacity_100C_air,viscosity_100C_air,sum(x*M),method)

    T=273.15+500 # K
    p=1e5 # Pa
    viscosity=mu(T)
    density=1e5/8.3145/T*M/1000
    heat_capacity=cp(T).squeeze()
    conductivity=conductivity_fun(heat_capacity,viscosity,M,method)

    viscosity_500=viscosity;conductivity_500=conductivity # store for comparison

    # air properties 500°C
    x=array([0,0,0,0,0,0.7812,0.2096,0.0092,0,0,0,0]) # N2,O2,Ar
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_500C_air=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_500C_air=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_500C_air=heat_capacity.dot(x)
    density_500C_air=1e5/8.3145/T*M.dot(x)/1000
    conductivity_500C_air_2=conductivity_fun(heat_capacity_500C_air,viscosity_500C_air,sum(x*M),method)

    if method.lower()=='eucken':
        # only print viscosity error table once
        print('error in viscosity estimation by pure monatomic gas assumption (Chapman-Enskog)')
        print('\tmu(100°C) table/(Pa s)\tmu(100°C) calc/(Pa s)\t% error(100°C)\t\tmu(500°C) table/(Pa s)\tvisc(500°C) calc/(Pa s)\t% error(500°C)')
        for i in range(len(M)):
            print('{:s}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}'.format(labels[i],
                                                                            visc_tab_100C[i],viscosity_100[i],(viscosity_100[i]/visc_tab_100C[i]-1)*100,
                                                                            visc_tab_500C[i],viscosity_500[i],(viscosity_500[i]/visc_tab_500C[i]-1)*100))
        print('{:s}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}'.format('air',
                                                                            visc_tab_100C_air,viscosity_100C_air,(viscosity_100C_air/visc_tab_100C_air-1)*100,
                                                                            visc_tab_500C_air,viscosity_500C_air,(viscosity_500C_air/visc_tab_500C_air-1)*100))

    print('\n')
    print('error in conductivity estimation by pure monatomic gas assumption and {:s} correlation'.format(method))
    print('\tk(100°C) table/(W/(Km))\tk(100°C) calc/(W/(Km))\t% error(100°C)\t\tk(500°C) table/(W/(Km))\tk(500°C) calc/(W/(Km))\t% error(500°C)')
    for i in range(len(M)):
        print('{:s}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}'.format(labels[i],
                                                                        conduct_tab_100C[i],conductivity_100[i],(conductivity_100[i]/conduct_tab_100C[i]-1)*100,
                                                                        conduct_tab_500C[i],conductivity_500[i],(conductivity_500[i]/conduct_tab_500C[i]-1)*100))
    print('{:s}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}'.format('air',
                                                                        conduct_tab_100C_air,conductivity_100C_air,(conductivity_100C_air/conduct_tab_100C_air-1)*100,
                                                                        conduct_tab_500C_air,conductivity_500C_air,(conductivity_500C_air/conduct_tab_500C_air-1)*100))
    print('{:s}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}\t{:16.16g}'.format('air',
                                                                        conduct_tab_100C_air,conductivity_100C_air_2,(conductivity_100C_air_2/conduct_tab_100C_air-1)*100,
                                                                        conduct_tab_500C_air,conductivity_500C_air_2,(conductivity_500C_air_2/conduct_tab_500C_air-1)*100))
    sse500=sqrt((conductivity_500/conduct_tab_500C-1).dot((conductivity_500/conduct_tab_500C-1)))
    sse100=sqrt((conductivity_100/conduct_tab_100C-1).dot((conductivity_100/conduct_tab_100C-1)))
    print(('SSE:\t'+'\t'*2+'{:'+str(16*3)+'.16g}'+'\t'*3+'{:'+str(16*3)+'.16g}').format(sse100,sse500))
```

    error in viscosity estimation by pure monatomic gas assumption (Chapman-Enskog)
    	mu(100°C) table/(Pa s)	mu(100°C) calc/(Pa s)	% error(100°C)		mu(500°C) table/(Pa s)	visc(500°C) calc/(Pa s)	% error(500°C)
    H2	        1.04e-05	1.02702531043997e-05	-1.247566303849001	        1.71e-05	1.655255653381486e-05	-3.20142377886049
    CH4	        1.35e-05	1.336496147017774e-05	-1.000285406090795	        2.34e-05	2.289723965628876e-05	-2.148548477398449
    CO2	        1.84e-05	1.843157288890713e-05	0.1715917875387829	        3.33e-05	3.264505253886455e-05	-1.966809192598951
    CO	        2.09e-05	2.076070611610776e-05	-0.6664779133599752	        3.45e-05	3.407368651020189e-05	-1.23569127477714
    H2O	        1.24e-05	1.321918149108586e-05	6.606302347466642	        2.86e-05	2.769240096431357e-05	-3.173423201700809
    N2	        2.11e-05	2.061962864551919e-05	-2.276641490430387	         3.5e-05	3.340860692534334e-05	-4.546837356161881
    O2	        2.47e-05	2.427683323183883e-05	-1.713225782028993	        4.14e-05	4.02902883427965e-05	-2.680462940105066
    Ar	         2.7e-05	2.682794796263007e-05	-0.6372297680367733	        4.65e-05	4.408212357838636e-05	-5.199734240029342
    He	        2.32e-05	2.330170909444836e-05	0.4384012691739825	        3.86e-05	3.737974558702046e-05	-3.16128086264128
    C6H6	         9.4e-06	9.488919167160711e-06	0.94594858681607	        1.93e-05	1.90423520273147e-05	-1.334963589043003
    C6H12	8.099999999999999e-06	8.692924624053638e-06	7.320057087081966	        1.54e-05	1.646119610923723e-05	6.89088382621581
    C6H14	8.099999999999999e-06	8.191795681180423e-06	1.133280014573135	        1.61e-05	1.635266780016617e-05	1.569365218423413
    air	        2.19e-05	2.144374724680555e-05	-2.08334590499748	       3.653e-05	3.493000582411623e-05	-4.379945731956647
    
    
    error in conductivity estimation by pure monatomic gas assumption and Eucken correlation
    	k(100°C) table/(W/(Km))	k(100°C) calc/(W/(Km))	% error(100°C)		k(500°C) table/(W/(Km))	k(500°C) calc/(W/(Km))	% error(500°C)
    H2	          0.2149	0.203382697859644	-5.359377450142389	          0.3678	0.3307648644270865	-10.06936801873669
    CH4	          0.0453	0.04123886384188439	-8.964980481491413	          0.1184	0.1031790327181253	-12.85554669077259
    CO2	          0.0225	0.0212210659729256	-5.684151231441781	          0.0536	0.04554189932489457	-15.03376991624148
    CO	          0.0302	0.02942703572270776	-2.559484361894826	          0.0534	0.05123273550965352	-4.058547734731254
    H2O	          0.0239	0.03264137934089893	36.5748089577361	          0.0669	0.07504522565433965	12.17522519333281
    N2	          0.0309	0.02917109290671209	-5.595168586692257	          0.0551	0.04968764919901746	-9.822778223198814
    O2	          0.0324	0.03056504853620499	-5.663430443811746	          0.0597	0.05535210002396424	-7.282914532723217
    Ar	          0.0209	0.02091187223034826	0.05680492989599273	          0.0362	0.03436117205824733	-5.079635198211796
    He	          0.1793	0.1816323648795678	1.300816999201215	          0.2983	0.2913679662744002	-2.323846371303995
    C6H6	          0.0169	0.01394504404976789	-17.48494645107753	          0.0618	0.04822014756211905	-21.97387125870704
    C6H12	          0.0203	0.015243303006264	-24.9098374075665	          0.0735	0.05558978353024756	-24.3676414554455
    C6H14	          0.0203	0.01723687480405856	-15.08928667951447	0.07530000000000001	0.05743311011606361	-23.72760940761805
    air	         0.03162	0.02938818069921645	-7.058252058139002	          0.0558	0.05071062740455404	-9.120739418361934
    air	         0.03162	0.02939024274938477	-7.051730710358095	          0.0558	0.05071804075332021	-9.107453847096403
    SSE:			                              0.5200631326616114			                              0.4977692836391809
    
    
    error in conductivity estimation by pure monatomic gas assumption and Eucken_mod correlation
    	k(100°C) table/(W/(Km))	k(100°C) calc/(W/(Km))	% error(100°C)		k(500°C) table/(W/(Km))	k(500°C) calc/(W/(Km))	% error(500°C)
    H2	          0.2149	0.2172301798654279	1.084308918300581	          0.3678	0.3540342535225912	-3.742726067810986
    CH4	          0.0453	0.04610111483102357	1.768465410648057	          0.1184	0.1219179549493422	2.971245734241701
    CO2	          0.0225	0.02383229916302765	5.921329613456239	          0.0536	0.0527127779542598	-1.655265010709339
    CO	          0.0302	0.03144593937130952	4.125627057316295	          0.0534	0.05548559404052496	3.905606817462459
    H2O	          0.0239	0.03575926137330548	49.62034047408148	          0.0669	0.08370986902208194	25.12685952478615
    N2	          0.0309	0.03115836556942978	0.8361345288989686	          0.0551	0.05368307065390265	-2.57155961179194
    O2	          0.0324	0.03277653322804318	1.162139592725864	          0.0597	0.06050256842035291	1.344335712483935
    Ar	          0.0209	0.02091187223034826	0.05680492989599273	          0.0362	0.03436117205824733	-5.079635198211796
    He	          0.1793	0.1816323648795678	1.300816999201215	          0.2983	0.2913679662744002	-2.323846371303995
    C6H6	          0.0169	0.01719368485797	1.737780224674546	          0.0618	0.06121479594979987	-0.9469321200649361
    C6H12	          0.0203	0.01908863144213823	-5.967332797348623	          0.0735	0.07142328569277905	-2.825461642477489
    C6H14	          0.0203	0.02180229736042627	7.400479608011201	0.07530000000000001	0.07391453821183336	-1.839922693448404
    air	         0.03162	0.03140315639283449	-0.6857799088093364	          0.0558	0.05490506120715654	-1.603832962085061
    air	         0.03162	0.03140442092343881	-0.6817807607880977	          0.0558	0.05490900482089199	-1.596765553957014
    SSE:			                              0.5114396746255642			                              0.2690868593184603
    
    
    error in conductivity estimation by pure monatomic gas assumption and nonpolar correlation
    	k(100°C) table/(W/(Km))	k(100°C) calc/(W/(Km))	% error(100°C)		k(500°C) table/(W/(Km))	k(500°C) calc/(W/(Km))	% error(500°C)
    H2	          0.2149	0.2098737050498552	-2.338899464934752	          0.3678	0.3416723905656043	-7.103754604240276
    CH4	          0.0453	0.04351804399304337	-3.933677719551065	          0.1184	0.1119629025140082	-5.436737741547148
    CO2	          0.0225	0.02244508153078594	-0.2440820853958248	          0.0536	0.04890324868240951	-8.762595741773293
    CO	          0.0302	0.03037339680798984	0.5741616158603913	          0.0534	0.0532262629459995	-0.3253502883904469
    H2O	          0.0239	0.03410288654358951	42.68990185602304	          0.0669	0.07910677723296881	18.24630378620151
    N2	          0.0309	0.03010262696736101	-2.580495251258852	          0.0551	0.05156050300599487	-6.423769499101873
    O2	          0.0324	0.03160168198550414	-2.463944489184733	          0.0597	0.05776638208477142	-3.238890980282372
    Ar	          0.0209	0.02091187223034826	0.05680492989597052	          0.0362	0.03436117205824733	-5.079635198211796
    He	          0.1793	0.1816323648795677	1.300816999201193	          0.2983	0.2913679662744001	-2.323846371304006
    C6H6	          0.0169	0.01546784442861263	-8.474293321818749	          0.0618	0.05431138899384442	-12.11749353746858
    C6H12	          0.0203	0.01704580071058004	-16.03053837152688	          0.0735	0.06301173766893418	-14.26974466811675
    C6H14	          0.0203	0.01937691662735592	-4.547208732236829	0.07530000000000001	0.06515877953595567	-13.46775626035105
    air	         0.03162	0.03033270055559991	-4.071155738140719	          0.0558	0.05267676824952395	-5.597189516982171
    air	         0.03162	0.0303343887684726	-4.065816671497158	          0.0558	0.05268255515999447	-5.586818709687325
    SSE:			                              0.4698569435587342			                               0.332360631216458



```python
import seaborn as sns
sns.set_context("talk", font_scale=0.8)
T=linspace(650,800)+273.15 # K
viscosity=mu(T)[:,:6+1]
density=101325/8.3145*outer(1/T,M[:6+1])/1000
heat_capacity=cp(T)[:,:6+1]
conductivity=(heat_capacity+5/4*R)*viscosity/M[:6+1]*1000

base_cm_1 = plt.get_cmap('viridis') # get colormap viridis
base_cm_2 = plt.get_cmap('inferno') # get colormap inferno
# join both selected color maps and split result into 18 parts
color_samples = array([
    *base_cm_1(linspace(0, 1, 8))[:, :3],
    *base_cm_2(linspace(0, 1, 10))[:, :3]
])
plt.rcParams['axes.prop_cycle'] = cycler(
    linestyle=['-','-.','--'])*cycler(
    fillstyle=['none','full'])*cycler(
    marker=['o', 'v','s'])+cycler(
    color=color_samples) # set color and marker cycle

fig,ax=plt.subplots(2,2,sharex=True)
fig.set_size_inches(10*16/9,10)

ax[0,0].plot(T-273.15,heat_capacity)
ylims=ax[0,0].get_ylim()
ax[0,0].vlines(1000-273.15,ylims[0],ylims[1],linestyle='--')
ax[0,0].text(1000-273.15,ylims[1],'1000K')
ax[0,0].set_xlabel(r'$T / ^\circ C$')
ax[0,0].set_ylabel(r'$Cp / (J/mol/K)$')
ax[0,0].legend(labels)

ax[0,1].plot(T-273.15,viscosity)
ylims=ax[0,1].get_ylim()
ax[0,1].vlines(1000-273.15,ylims[0],ylims[1],linestyle='--')
ax[0,1].text(1000-273.15,ylims[1],'1000K')
ax[0,1].set_xlabel(r'$T / ^\circ C$')
ax[0,1].set_ylabel(r'$\mu / (Pa s)$')
ax[0,1].legend(labels)

ax[1,0].semilogy(T-273.15,density)
ylims=ax[1,0].get_ylim()
ax[1,0].vlines(1000-273.15,ylims[0],ylims[1],linestyle='--')
ax[1,0].text(1000-273.15,ylims[1],'1000K')
ax[1,0].set_xlabel(r'$T / ^\circ C$')
ax[1,0].set_ylabel(r'$\rho / (kg/m^3)$')
ax[1,0].legend(labels)

ax[1,1].plot(T-273.15,conductivity)
ylims=ax[1,1].get_ylim()
ylims=ax[1,1].set_ylim(ylims[0],1)
ax[1,1].vlines(1000-273.15,ylims[0],ylims[1],linestyle='--')
ax[1,1].text(1000-273.15,ylims[1],'1000K')
ax[1,1].set_yscale('log')
ax[1,1].set_xlabel(r'$T / ^\circ C$')
ax[1,1].set_ylabel(r'$\lambda / (W/m/K)$')
ax[1,1].legend(labels[:6+1])
#ax[1,1].yaxis.set_major_locator(plt.LogLocator(base=10, subs=[0.1,1]))

print('conductivity values range between: {:0.5g} and {:0.5g} W/m/K'.format(conductivity.min(),conductivity.max()))
```

# Effective heat conductivity of metal foam

Equivalence to VDI-Wärmeatlas nomenclature

$\lambda_p$: heat conductivity of particles

$\lambda_f$: heat conductivity of fluid filling void spaces

$\psi$: Porosity of the bed

The analytical solution by Maxwell (1873) is valid for diluted suspensions or emulsions ($\psi$ is low), where particles do not interact thermically [4]. 

$\lambda_{bed}/\lambda_{f}=\frac{1+2\frac{(1-\psi)\left(\frac{\lambda_p}{\lambda_f}-1\right)}{\left(\frac{\lambda_p}{\lambda_f}+2\right)}}{1-\frac{(1-\psi)\left(\frac{\lambda_p}{\lambda_f}-1\right)}{\left(\frac{\lambda_p}{\lambda_f}+2\right)}}=\frac{\lambda_p+2\lambda_f+2(1-\psi)(\lambda_p-\lambda_f)}{\lambda_p+2\lambda_f-(1-\psi)(\lambda_p-\lambda_f)}$

For beds, values are acceptable if these conditions are fulfilled:

1. Particles do not influence one another ==>> Diluted bed: $\psi<=0.40$
2. Moderate particle-to-fluid conductivity ratio: $\lambda_p/\lambda_f<20$.

Since $\lambda_{Ni,1000K}=72\frac{W}{mK}$ [5], $\lambda_{H_2,1000K}=0.40\frac{W}{mK}$ and $\lambda_{N_2,1000K}=0.061\frac{W}{mK}$, then $\rightarrow \lambda_{fuel}/\lambda_{Ni}\approx 0.0055 \text{ to } 0.00084$ and it is best to choose $\lambda_p=\lambda_{fuel}$ and $\lambda_f=\lambda_{Ni}$. In this case the porosity of the fuel channel is $\epsilon_{FC}=(1-\psi)$. The "diluted bed" condition is also fulfilled, since at low values of $\psi=\frac{V_f}{V_f+V_p}=\frac{V_{Ni}}{V_{Ni}+V_{fuel}}$, the value of $\epsilon_{FC}=(1-\psi)=\frac{V_{fuel}}{V_{Ni}+V_{fuel}}\approx 93\%$ is high.

$\lambda_{bed}/\lambda_{Ni}=\frac{\lambda_{fuel}+2\lambda_{Ni}+2\epsilon_{FC}(\lambda_{fuel}-\lambda_{Ni})}{\lambda_{fuel}+2\lambda_{Ni}-\epsilon_{FC}(\lambda_{fuel}-\lambda_{Ni})}$


```python
import seaborn as sns
sns.set_context("talk", font_scale=0.8)
efc=linspace(0.6,1,30) # porosity of the fuel channel
k_Ni=72 # thermal conductivity of nickel at 1000K in W/m/K
fuel_labels=['H_2','CO_2']
k_fuel_values=[0.4,0.057] # thermal conductivity of fuel at 1000K in W/m/K

# plt.rcParams['axes.prop_cycle'] = cycler(
#     linestyle=['--','-.',':','-'])*cycler(color=['b', 'orange', 'g', 'r']) # set color and marker cycle

fig,ax=plt.subplots(2,2,sharex=True)
fig.set_size_inches(10*16/9,10)

for i in range(len(k_fuel_values)):
    k_fuel=k_fuel_values[i]
    k_Ni_foam=k_Ni*(k_fuel+2*k_Ni+2*efc*(k_fuel-k_Ni))/(k_fuel+2*k_Ni-efc*(k_fuel-k_Ni)) # Maxwell 1873 approach
    ax[0,0].plot(efc,k_Ni_foam,label='$'+r'\lambda_{fuel}=\lambda_{'+fuel_labels[i]+'}='+'{:g}~W/(m\cdot K)'.format(k_fuel)+'$')
# determine required k_fuel for a foam conductivity of 8.46W/m/K at porosity 0.867 by Maxwell approach
min_k_Ni_foam=min(k_Ni_foam)
k_fuel=72*((2-2*0.867-8.46/72*(2+0.867))/(-1-2*0.867+8.46/72*(1-0.867)))
k_Ni_foam=k_Ni*(k_fuel+2*k_Ni+2*0.867*(k_fuel-k_Ni))/(k_fuel+2*k_Ni-0.867*(k_fuel-k_Ni)) # Maxwell 1873 approach
ax[0,0].plot(0.867,k_Ni_foam,'o',color='black',fillstyle='full')
ax[0,0].plot([min(efc),0.867,0.867],[k_Ni_foam,k_Ni_foam,min_k_Ni_foam],'--',color='black')
ax[0,0].text(0.867*(1+.01),min_k_Ni_foam,r'$\epsilon_{FC}='+'{:g}'.format(0.867)+'$')
ax[0,0].text(0.867*(1+.01),k_Ni_foam,r'$\lambda_{Fuel}='+'{:g}'.format(k_fuel)+r'~\frac{W}{m\cdot K}'+'$')
ax[0,0].text(min(efc),k_Ni_foam,r'$\lambda_{Ni-foam}='+'{:g}'.format(k_Ni_foam)+r'~\frac{W}{m\cdot K}'+'$',horizontalalignment='left',verticalalignment='bottom')
k_fuel=k_fuel_values[0]
k_Ni_foam=k_Ni*(k_fuel+2*k_Ni+2*0.867*(k_fuel-k_Ni))/(k_fuel+2*k_Ni-0.867*(k_fuel-k_Ni)) # Maxwell 1873 approach
ax[0,0].plot([min(efc),0.867,0.867],[k_Ni_foam,k_Ni_foam,min_k_Ni_foam],'--')
ax[0,0].text(min(efc),k_Ni_foam,r'$\lambda_{Ni-foam}='+'{:g}'.format(k_Ni_foam)+r'~\frac{W}{m\cdot K}'+'$',horizontalalignment='left',verticalalignment='top')
#ax[0,0].set_xlabel(r'$\epsilon_{FC}$ / -')
ax[0,0].set_ylabel(r'$\lambda_{Ni-foam}$ / (W/m/K)')
ax[0,0].set_title(r'Maxwell approach $\lambda_{Ni-foam}=\lambda_{Ni}\frac{\lambda_{fuel}+2\lambda_{Ni}+2\epsilon_{FC}(\lambda_{fuel}-\lambda_{Ni})}{\lambda_{fuel}+2\lambda_{Ni}-\epsilon_{FC}(\lambda_{fuel}-\lambda_{Ni})}$',pad=20)
ax[0,0].legend()

for i in range(len(k_fuel_values)):
    k_fuel=k_fuel_values[i]
    k_Ni_foam=k_Ni/((1-efc)+efc*k_Ni/k_fuel)
    ax[1,0].plot(efc,k_Ni_foam,label='$'+fuel_labels[i]+'$')
ax[1,0].set_xlabel(r'$\epsilon_{FC}$ / -')
ax[1,0].set_ylabel(r'$\lambda_{Ni-foam}$ / (W/m/K)')
ax[1,0].set_title(r'Series approach $\lambda_{Ni-foam,I}=\left(\frac{1-\epsilon_{FC}}{\lambda_{Ni}}+\frac{\epsilon_{FC}}{\lambda_{fuel}}\right)^{-1}$')
ax[1,0].legend()

for i in range(len(k_fuel_values)):
    k_fuel=k_fuel_values[i]
    k_Ni_foam=(1-efc)*k_Ni+efc*k_fuel
    ax[1,1].plot(efc,k_Ni_foam,label='$'+fuel_labels[i]+'$')
ax[1,1].set_xlabel(r'$\epsilon_{FC}$ / -')
ax[1,1].set_ylabel(r'$\lambda_{Ni-foam}$ / (W/m/K)')
ax[1,1].set_title(r'Parallel approach $\lambda_{Ni-foam,II}=(1-\epsilon_{FC})\lambda_{Ni}+\epsilon_{FC}\lambda_{fuel}$')
ax[1,1].legend()

for a in [0.05,0.2,1-.35]: # series to parallel proportion
    for i in range(len(k_fuel_values)):
        k_fuel=k_fuel_values[i]
        k_Ni_foam_series=k_Ni/((1-efc)+efc*k_Ni/k_fuel)
        k_Ni_foam_parallel=(1-efc)*k_Ni+efc*k_fuel
        k_Ni_foam=(a/k_Ni_foam_series+(1-a)/k_Ni_foam_parallel)**-1
        ax[0,1].plot(efc,k_Ni_foam,label='a={:g}; '.format(a)+'$'+fuel_labels[i]+'$')
#ax[0,1].set_xlabel(r'$\epsilon_{FC}$ / -')
ax[0,1].set_ylabel(r'$\lambda_{Ni-foam}$ / (W/m/K)')
ax[0,1].set_title(r'$\lambda_{Ni-foam,combined}=\left(\frac{a}{\lambda_{Ni-foam,I}}+\frac{1-a}{\lambda_{Ni-foam,II}}\right)^{-1}$',pad=20)
ax[0,1].legend(loc='upper right')

plt.savefig('summary_series_parallel_approach.png',dpi=80)
```

Some values at 1000K

$\begin{align}
\lambda_{Ni,1000K}&=72\frac{W}{K m}\\
\lambda_{Ni-foam,1000K}&=8.46\frac{W}{K m}\\
\epsilon_{FC}&=0.867\\
\lambda_{Fuel,1000K}&=? \leftarrow \text{solve from } \lambda_{Ni-foam}=\lambda_{Ni}\frac{\lambda_{fuel}+2\lambda_{Ni}+2\epsilon_{FC}(\lambda_{fuel}-\lambda_{Ni})}{\lambda_{fuel}+2\lambda_{Ni}-\epsilon_{FC}(\lambda_{fuel}-\lambda_{Ni})}
\end{align}$



```python
efc=0.867 # porosity of the fuel channel
k_Ni=72 # thermal conductivity of nickel at 1000K in W/m/K
k_fuel=72*((2-2*efc-8.46/72*(2+efc))/(-1-2*efc+8.46/72*(1-efc)))
print('lambda_fuel={:0.4g} W/(m K)'.format(k_fuel))
print('confirm lambda_Ni-foam={:0.4g} W/(m K)'.format(k_Ni*(k_fuel+2*k_Ni+2*efc*(k_fuel-k_Ni))/(k_fuel+2*k_Ni-efc*(k_fuel-k_Ni))))
```

    lambda_fuel=1.877 W/(m K)
    confirm lambda_Ni-foam=8.46 W/(m K)


Assuming combined series and parallel approach $\lambda_{Ni-foam,combined}=\left(\frac{a}{\lambda_{Ni-foam,I}}+\frac{1-a}{\lambda_{Ni-foam,II}}\right)^{-1}$, 
Which value of a is needed to obtain $\lambda_{eff}=8.46 \frac{W}{m K}?$ with different fuel compositions.


```python
efc=0.867 # porosity of the fuel channel
T=1000 # K
p=1e5 # Pa
viscosity=mu(T)
density=1e5/8.3145/T*M/1000
heat_capacity=cp(T).squeeze()
conductivity=conductivity_fun(heat_capacity,viscosity,M,'Eucken_mod')

viscosity_100=viscosity;conductivity_100=conductivity # store for comparison

# 1000K
# labels=['H2','CH4','CO2','CO','H2O','N2','O2','Ar','He','C6H6','C6H12','C6H14'] # order of components
x_compositions=array([
    [0.4,0,0,0,0,0.6,0,0,0,0,0,0], # Santhanam 2017 SOFC 40%H2 in N2
    [0.1,0,0,0,0.9,0,0,0,0,0,0,0], # Santhanam 2017 SOEC 10%H2, 90%H2O Santhanam 2017
    [0.448,0.069,0.117,0.054,0.312,0,0,0,0,0,0,0], # Mottaghizadeh 2017 SOFC
    [0.078,0,0.185,0.025,0.712,0,0,0,0,0,0,0] # Mottaghizadeh 2017 SOEC
]) # consider different fuel/reagent compositions

print('porosity of the fuel channel: {:0.4g}'.format(efc))
print('Conductivity of Nickel at 1000K: {:0.4g} W/(m K)'.format(k_Ni))
print('Conductivity of fuel at 1000K and corresponding Ni-foam value assuming series/parallel approach')
for x in x_compositions: # loop each line of the compositions
    # calculate conductivity of gas mixture
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_mix=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_mix=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_mix=heat_capacity.dot(x)
    density_mix=1e5/8.3145/T*M.dot(x)/1000
    
    k_fuel=conductivity_mix
    
    k_maxwell=k_Ni*(k_fuel+2*k_Ni+2*efc*(k_fuel-k_Ni))/(k_fuel+2*k_Ni-efc*(k_fuel-k_Ni))
    
    # calculate conductivity of foam by series and parallel approach
    k_Ni_foam_series=k_Ni/((1-efc)+efc*k_Ni/k_fuel)
    k_Ni_foam_parallel=(1-efc)*k_Ni+efc*k_fuel
    
    # calculate a parameter to meet 8.46 W/m/K
    # k_Ni_foam=(a/k_Ni_foam_series+(1-a)/k_Ni_foam_parallel)**-1
    k_eff=8.46 # W/m/K assumption
    a=(1/k_eff-1/k_Ni_foam_parallel)/(1/k_Ni_foam_series+1/k_Ni_foam_parallel)
    
    gas_name=':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0])
    print(('{:30s}\tk_gas={:5.4g} W/(m K) \t-->>k_series={:5.4g} W/(m K), k_parallel={:5.4g} W/(m K). k_eff={:5.4g} W/(m K) <<--a={:5.4g} Maxwell: {:5.4g} W/(m K)').format(gas_name,conductivity_mix,k_Ni_foam_series,k_Ni_foam_parallel,k_eff,a,k_maxwell))
```

    porosity of the fuel channel: 0.867
    Conductivity of Nickel at 1000K: 72 W/(m K)
    Conductivity of fuel at 1000K and corresponding Ni-foam value assuming series/parallel approach
    H2:N2 40:60                   	k_gas=0.1668 W/(m K) 	-->>k_series=0.1924 W/(m K), k_parallel=9.721 W/(m K). k_eff= 8.46 W/(m K) <<--a=0.002892 Maxwell: 6.839 W/(m K)
    H2:H2O 10:90                  	k_gas=0.1371 W/(m K) 	-->>k_series=0.158 W/(m K), k_parallel=9.695 W/(m K). k_eff= 8.46 W/(m K) <<--a=0.002341 Maxwell:  6.81 W/(m K)
    H2:CH4:CO2:CO:H2O 45:7:12:5:31	k_gas=0.2129 W/(m K) 	-->>k_series=0.2455 W/(m K), k_parallel=9.761 W/(m K). k_eff= 8.46 W/(m K) <<--a=0.003772 Maxwell: 6.882 W/(m K)
    H2:CO2:CO:H2O 8:18:2:71       	k_gas=0.1188 W/(m K) 	-->>k_series=0.137 W/(m K), k_parallel=9.679 W/(m K). k_eff= 8.46 W/(m K) <<--a=0.002011 Maxwell: 6.793 W/(m K)



```python
T=273.15+726.5#(52+25)/2 # K
p=101325 # Pa
# labels=['H2','CH4','CO2','CO','H2O','N2','O2','Ar','He','C6H6','C6H12','C6H14'] # order of components
x_compositions=array([
    [0.4,0,0,0,0,0.6,0,0,0,0,0,0], # Santhanam 2017 SOFC 40%H2 in N2
    [0.1,0,0,0,0.9,0,0,0,0,0,0,0], # Santhanam 2017 SOEC 10%H2, 90%H2O Santhanam 2017
    [0.448,0.069,0.117,0.054,0.312,0,0,0,0,0,0,0], # Mottaghizadeh 2017 SOFC
    [0.078,0,0.185,0.025,0.712,0,0,0,0,0,0,0], # Mottaghizadeh 2017 SOEC
    [0,0,0,0,0,0.7812,0.2096,0.0092,0,0,0,0], # air
    [0,0,0,0,0,0.79,0.21,0,0,0,0,0], # air
    [0.25,0,0,0,0.15,0.6,0,0,0,0,0,0], # additional for validation - Marius
    [40/240,0,0,0,0,200/240,0,0,0,0,0,0], # info. H2:N2 40:200
    [50/250,0,0,0,0,200/250,0,0,0,0,0,0], # info. H2:N2 50:200
    [0.4918886373731  ,0.07979978711003 ,0.07991435965461 ,0.06158275027688 ,0.2868144655854  ,0.0,0,0,0,0,0,0              ] # switch with methane
]) # consider different fuel/reagent compositions
print('Cp,rho,k of representative gases at {:0.4g}K'.format(T))
for x in x_compositions: # loop each line of the compositions
    # calculate conductivity of gas mixture
    viscosity=mu(T)
    density=p/8.3145/T*M/1000
    heat_capacity=cp(T).squeeze()
    conductivity=(heat_capacity+5/4*R)*viscosity/M*1000
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_mix=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_mix=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_mix=heat_capacity.dot(x)/M.dot(x)*1000
    density_mix=p/8.3145/T*M.dot(x)/1000
    heat_capacity
    gas_name=':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0])
    print(('{:30s}\tk={:5.4g} W/(m K) \trho={:5.4g} kg/m^3\tCp={:5.4g} J/(kg K)\t eta={:5.4g} Pa s').format(gas_name,conductivity_mix,density_mix,heat_capacity_mix,viscosity_mix))
```

    Cp,rho,k of representative gases at 999.6K
    H2:N2 40:60                   	k=0.1545 W/(m K) 	rho=0.2146 kg/m^3	Cp= 1804 J/(kg K)	 eta=3.843e-05 Pa s
    H2:H2O 10:90                  	k=0.1227 W/(m K) 	rho=0.1999 kg/m^3	Cp= 2449 J/(kg K)	 eta=3.562e-05 Pa s
    H2:CH4:CO2:CO:H2O 45:7:12:5:31	k=0.1933 W/(m K) 	rho=0.174 kg/m^3	Cp= 2771 J/(kg K)	 eta=3.545e-05 Pa s
    H2:CO2:CO:H2O 8:18:2:71       	k=0.1057 W/(m K) 	rho=0.2659 kg/m^3	Cp= 1954 J/(kg K)	 eta=3.734e-05 Pa s
    N2:O2:Ar 78:21:1              	k=0.06196 W/(m K) 	rho=0.3529 kg/m^3	Cp= 1142 J/(kg K)	 eta=4.129e-05 Pa s
    N2:O2 79:21                   	k=0.06215 W/(m K) 	rho=0.3516 kg/m^3	Cp= 1150 J/(kg K)	 eta=4.116e-05 Pa s
    H2:H2O:N2 25:15:60            	k=0.1224 W/(m K) 	rho=0.2438 kg/m^3	Cp= 1670 J/(kg K)	 eta=3.849e-05 Pa s
    H2:N2 17:83                   	k=0.09541 W/(m K) 	rho=0.2885 kg/m^3	Cp= 1365 J/(kg K)	 eta=3.928e-05 Pa s
    H2:N2 20:80                   	k=0.103 W/(m K) 	rho=0.278 kg/m^3	Cp= 1414 J/(kg K)	 eta=3.921e-05 Pa s
    H2:CH4:CO2:CO:H2O 49:8:8:6:29 	k=0.2085 W/(m K) 	rho=0.1544 kg/m^3	Cp= 3069 J/(kg K)	 eta=3.462e-05 Pa s



```python
T=273.15+702.8#(52+25)/2 # K
p=101325 # Pa
# labels=['H2','CH4','CO2','CO','H2O','N2','O2','Ar','He','C6H6','C6H12','C6H14'] # order of components
x_compositions=array([
    [0.4918886373731  ,0.07979978711003 ,0.07991435965461 ,0.06158275027688 ,0.2868144655854  ,0.0,0,0,0,0,0,0              ] # switch with methane
]) # consider different fuel/reagent compositions
print('Cp,rho,k,Dim of gas mixture at {:0.4g} K and {:7g} Pa'.format(T,p))
for x in x_compositions: # loop each line of the compositions
    # calculate conductivity of gas mixture
    viscosity=mu(T)
    density=p/8.3145/T*M/1000
    heat_capacity=cp(T).squeeze()
    conductivity=(heat_capacity+5/4*R)*viscosity/M*1000
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_mix=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_mix=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_mix=heat_capacity.dot(x)/M.dot(x)*1000
    density_mix=p/8.3145/T*M.dot(x)/1000
    diffusivity_im=(1-x)/(x/D_ij(T,p)-diag(diag(x/D_ij(T,p)))).sum(axis=1)
    gas_name=':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0])
    print(('{:30s}\tk={:5.4g} W/(m K) \trho={:5.4g} kg/m^3\tCp={:5.4g} J/(kg K)\t eta={:5.4g} Pa s').format(gas_name,conductivity_mix,density_mix,heat_capacity_mix,viscosity_mix))
    print('Diffusivities of each component in the mixture Dim/(m^2/s)')
    print(('{:10s}'*len(labels)).format(*labels))
    print(('{:<10.4g}'*len(diffusivity_im)).format(*diffusivity_im))
```

    Cp,rho,k,Dim of gas mixture at 975.9 K and  101325 Pa
    H2:CH4:CO2:CO:H2O 49:8:8:6:29 	k=0.2042 W/(m K) 	rho=0.1581 kg/m^3	Cp= 3051 J/(kg K)	 eta=3.397e-05 Pa s
    Diffusivities of each component in the mixture Dim/(m^2/s)
    H2        CH4       CO2       CO        H2O       N2        O2        Ar        He        C6H6      C6H12     C6H14     
    0.0005845 0.0002757 0.000224  0.0002636 0.0003505 0.0002503 0.0002593 0.000248  0.0007381 0.0001185 0.000102  0.0001028 



```python
print('Binary diffusion coefficients at {:0.4g} K and {:7g} Pa'.format(T,p))
print('D_ij/(m^2/s)')
print('{:17s}'.format('i (down) / j->')+('{:<10s}'*len(labels)).format(*labels))
print('\n'.join(['{:17s}'.format(labels[i])+''.join(['{:<10.3g}'.format(cell) for cell in row]) for i,row in enumerate(D_ij(T,p))]))
```

    Binary diffusion coefficients at 975.9 K and  101325 Pa
    D_ij/(m^2/s)
    i (down) / j->   H2        CH4       CO2       CO        H2O       N2        O2        Ar        He        C6H6      C6H12     C6H14     
    H2               0.00103   0.000532  0.000474  0.000551  0.000654  0.000544  0.000581  0.00057   0.00114   0.000301  0.000255  0.000262  
    CH4              0.000532  0.000178  0.000136  0.000167  0.000197  0.000166  0.000172  0.000164  0.000496  8.45e-05  7.31e-05  7.4e-05   
    CO2              0.000474  0.000136  9.23e-05  0.00012   0.000146  0.000119  0.000121  0.000113  0.000429  5.4e-05   4.67e-05  4.7e-05   
    CO               0.000551  0.000167  0.00012   0.000151  0.000188  0.00015   0.000154  0.000145  0.000504  7.26e-05  6.23e-05  6.32e-05  
    H2O              0.000654  0.000197  0.000146  0.000188  0.000205  0.000187  0.000193  0.000185  0.000634  8.12e-05  7.02e-05  7e-05     
    N2               0.000544  0.000166  0.000119  0.00015   0.000187  0.000149  0.000153  0.000144  0.000496  7.27e-05  6.24e-05  6.34e-05  
    O2               0.000581  0.000172  0.000121  0.000154  0.000193  0.000153  0.000157  0.000147  0.000531  7.15e-05  6.11e-05  6.21e-05  
    Ar               0.00057   0.000164  0.000113  0.000145  0.000185  0.000144  0.000147  0.000137  0.000518  6.61e-05  5.64e-05  5.73e-05  
    He               0.00114   0.000496  0.000429  0.000504  0.000634  0.000496  0.000531  0.000518  0.00118   0.000268  0.000225  0.000232  
    C6H6             0.000301  8.45e-05  5.4e-05   7.26e-05  8.12e-05  7.27e-05  7.15e-05  6.61e-05  0.000268  3.12e-05  2.77e-05  2.74e-05  
    C6H12            0.000255  7.31e-05  4.67e-05  6.23e-05  7.02e-05  6.24e-05  6.11e-05  5.64e-05  0.000225  2.77e-05  2.47e-05  2.45e-05  
    C6H14            0.000262  7.4e-05   4.7e-05   6.32e-05  7e-05     6.34e-05  6.21e-05  5.73e-05  0.000232  2.74e-05  2.45e-05  2.42e-05  



```python
print('x_i/D_ij(T,p) for i!=j at {:0.4g} K and {:7g} Pa'.format(T,p))
print('x_i/D_ij(T,p) for i!=j /(s/m^2)')
print('{:17s}'.format('i (down) / j->')+('{:<10s}'*len(labels)).format(*labels))
print('\n'.join(['{:17s}'.format(labels[i])+''.join(['{:<10.3g}'.format(cell) for cell in row]) for i,row in enumerate(x/D_ij(T,p)-diag(diag(x/D_ij(T,p))))]))
```

    x_i/D_ij(T,p) for i!=j at 975.9 K and  101325 Pa
    x_i/D_ij(T,p) for i!=j /(s/m^2)
    i (down) / j->   H2        CH4       CO2       CO        H2O       N2        O2        Ar        He        C6H6      C6H12     C6H14     
    H2               0         150       169       112       439       0         0         0         0         0         0         0         
    CH4              925       0         587       368       1.46e+03  0         0         0         0         0         0         0         
    CO2              1.04e+03  586       0         513       1.97e+03  0         0         0         0         0         0         0         
    CO               893       477       665       0         1.52e+03  0         0         0         0         0         0         0         
    H2O              753       406       549       327       0         0         0         0         0         0         0         0         
    N2               905       480       669       410       1.53e+03  0         0         0         0         0         0         0         
    O2               847       465       662       400       1.48e+03  0         0         0         0         0         0         0         
    Ar               863       486       709       424       1.55e+03  0         0         0         0         0         0         0         
    He               433       161       186       122       452       0         0         0         0         0         0         0         
    C6H6             1.64e+03  945       1.48e+03  849       3.53e+03  0         0         0         0         0         0         0         
    C6H12            1.93e+03  1.09e+03  1.71e+03  989       4.08e+03  0         0         0         0         0         0         0         
    C6H14            1.88e+03  1.08e+03  1.7e+03   974       4.09e+03  0         0         0         0         0         0         0         


Compare with Cantera data


```python
import cantera as ct
T=273.15+702.8#(52+25)/2 # K
p=101325 # Pa
gas=ct.Solution('gri30.yaml')
composition={labels[i]:x[i] for i in range(6) }
gas.TPX = T, 101325, composition
gas()

x_ax=arange(len(labels[:6]))
width=0.35

sns.set_context("talk", font_scale=0.8)
base_cm_1 = plt.get_cmap('viridis') # get colormap viridis
base_cm_2 = plt.get_cmap('inferno') # get colormap inferno
base_cm_1 = plt.get_cmap('tab10') # get colormap inferno
# join both selected color maps and split result into 18 parts
color_samples = array([
    *base_cm_1(linspace(0, 1, 8))[:, :3],
    *base_cm_2(linspace(0, 1, 10))[:, :3]
])
color_samples = array([
    *base_cm_1(linspace(0, 1, 10))[:, :3],
    *base_cm_1(linspace(0, 1, 8))[:, :3],
    #*base_cm_1(linspace(0, 1, 6))[:, :3]
])
plt.rcParams['axes.prop_cycle'] = cycler(
    linestyle=['-','-.','--'])*cycler(
    fillstyle=['none','full'])*cycler(
    marker=['o', 'v','s'])+cycler(
    color=color_samples) # set color and marker cycle

plt.bar(x_ax-width/2,gas[labels[:6]].mix_diff_coeffs_mole,width,label='cantera')#,color=plt.rcParams['axes.prop_cycle'].by_key()['color'][2])
plt.bar(x_ax+width/2,diffusivity_im[:6],width,label="Blanc's law [2,19]")#,color=plt.rcParams['axes.prop_cycle'].by_key()['color'][-3]);
plt.xticks(x_ax)
plt.gca().set_xticklabels(labels[:6])
plt.legend(loc='best')
plt.ylabel('$D_{im}$ / (m$^2$/s)')
plt.xlabel('component i in mixture');
```


```python
from numpy import finfo
```


```python
T=1000 # K
p=1e5 # Pa
t=linspace(650,850,20)+273.15

fig,ax_list=plt.subplots(ncols=2,nrows=1,figsize=(5*2,5))

h=dict({labels[i]:0 for i in range(6)})
s=dict({labels[i]:0 for i in range(6)})
g=dict({labels[i]:0 for i in range(6)})

delta_r_h=zeros(len(t))
delta_r_g=zeros(len(t))
delta_r_s=zeros(len(t))

print('H,S,G at {:0.4g} K and {:0.4g} bar'.format(T,p/1e5))

for component in labels[:6]:
    gas.TPX=T,p,component+':1'
    h[component]=gas.enthalpy_mole/1000 # J/mol
    s[component]=gas.entropy_mole/1000 # J/(mol K)
    g[component]=gas.gibbs_mole/1000 # J/mol
    delta_r_h[0]=-h['CH4']-h['H2O']+3*h['H2']+h['CO']
    delta_r_g[0]=-g['CH4']-g['H2O']+3*g['H2']+g['CO']
    delta_r_s[0]=-s['CH4']-s['H2O']+3*s['H2']+s['CO']
print('\n'.join(['{:7s}= {:<+10.4g} kJ/mol   {:7s}= {:<+10.4g} J/(mol K)   {:7s}= {:<+10.4g} kJ/mol'.format(
    'H('+labels[i]+')',h[labels[i]]/1000,
    'S('+labels[i]+')',s[labels[i]],
    'G('+labels[i]+')',g[labels[i]]/1000
) for i in range(6)]))

print('\n{:7s}= {:<+10.5g} kJ/mol   {:7s}= {:<+10.5g} J/(mol K)   {:7s}= {:<+10.5g} kJ/mol'.format(
    'delta_R_H('+'{:0.4g}'.format(T)+' K)',delta_r_h[0]/1000,
    'delta_R_S('+'{:0.4g}'.format(T)+' K)',delta_r_s[0],
    'delta_R_G('+'{:0.4g}'.format(T)+' K)',delta_r_g[0]/1000)
)


for i,T in enumerate(t):
    for component in labels[:6]:
        gas.TPX=T,p,component+':1'
        h[component]=gas.enthalpy_mole/1000 # J/mol
        s[component]=gas.entropy_mole/1000 # J/(mol K)
        g[component]=gas.gibbs_mole/1000 # J/mol
    
    delta_r_h[i]=-h['CH4']-h['H2O']+3*h['H2']+h['CO']
    delta_r_g[i]=-g['CH4']-g['H2O']+3*g['H2']+g['CO']
    delta_r_s[i]=-s['CH4']-s['H2O']+3*s['H2']+s['CO']

ax_list[0].plot(t,delta_r_g,label=r'$\Delta G(T)$ (cantera)')
ax_list[0].plot(t,+224.990e3-252.46*t,label=r'$(224990-252.46\cdot\frac{T}{K})\frac{J}{mol}$')
ax_list[0].set_ylabel(r'$\Delta G(T)$ / J mol$^{-1}$')
ax_list[0].set_xlabel(r'$T$ / K')


ax_list[1].plot(1/t,-delta_r_g/8.3145*1/t,label=r'$\left(\frac{-\Delta G(T)}{RT}\right)$ (cantera)')
ax_list[1].plot(1/t,-224.990e3/8.3145*1/t+252.46/8.3145,label=r'$-\frac{224990\frac{J}{mol}}{R}\cdot\frac{1}{T/K}+\frac{252.46\frac{J}{mol K}}{R}$')
ax_list[1].set_ylabel(r'$ln(K_p(T))$ / -')
ax_list[1].set_xlabel(r'$\frac{1}{T}$ / K$^{-1}$')
secondary_x=ax_list[1].secondary_xaxis('top',functions=(lambda x:1/(x+finfo('float').eps)-273.15,lambda x:1/(x+273.15+finfo('float').eps)))
secondary_x.set_xlabel('$T$ / °C')

for k in range(len(ax_list)):
    ax_list[k].legend(loc='best')
plt.tight_layout();
```

Parameter a can be varied. Which value is optimal? In ref. [6], a best fit for experimental data of high porosity metal foams (Al-water, Al-air) was found to be a=1-A=0.65. Value 0.2 in the Wärmeatlas is arbitrary.


```python
efc=0.867 # porosity of the fuel channel
A=0.35
print('porosity of the fuel channel: {:0.4g}'.format(efc))
print('Conductivity of Nickel at 1000K: {:0.4g} W/(m K)'.format(k_Ni))
print('Conductivity of fuel at 1000K and corresponding Ni-foam value assuming series/parallel approach with A={:0.4g}'.format(A))
for x in x_compositions: # loop each line of the compositions
    # calculate conductivity of gas mixture
    viscosity=mu(T)
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_mix=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_mix=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_mix=heat_capacity.dot(x)
    density_mix=1e5/8.3145/T*M.dot(x)/1000
    
    k_fuel=conductivity_mix
    
    # calculate conductivity of foam by series and parallel approach
    k_Ni_foam_series=k_Ni/((1-efc)+efc*k_Ni/k_fuel)
    k_Ni_foam_parallel=(1-efc)*k_Ni+efc*k_fuel
    
    # calculate parameter a based on A
    a=(k_Ni_foam_parallel/(A*k_Ni_foam_parallel+(1-A)*k_Ni_foam_series)-1)*k_Ni_foam_series/(k_Ni_foam_parallel-k_Ni_foam_series)
    
    # calculate k_eff based on a value
    k_eff=(a/k_Ni_foam_series+(1-a)/k_Ni_foam_parallel)**-1
    
    gas_name=':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0])
    print(('{:30s}\tk_gas={:5.4g} W/(m K) \t-->>k_series={:5.4g} W/(m K), k_parallel={:5.4g} W/(m K), a={:5.4g}. k_eff={:5.4g} W/(m K)').format(gas_name,conductivity_mix,k_Ni_foam_series,k_Ni_foam_parallel,a,k_eff))
```

    porosity of the fuel channel: 0.867
    Conductivity of Nickel at 1000K: 72 W/(m K)
    Conductivity of fuel at 1000K and corresponding Ni-foam value assuming series/parallel approach with A=0.35
    H2:CH4:CO2:CO:H2O 49:8:8:6:29 	k_gas=0.2051 W/(m K) 	-->>k_series=0.2365 W/(m K), k_parallel=9.754 W/(m K), a=0.04309. k_eff=3.568 W/(m K)


For some cells, the fuel GDL is Ni-foam with 93% porosity and the air GDL is Crofer 22 APU foam with 93% porosity.


```python
A=0.35 # 
efc=0.93
eac=0.93
k_Ni=72 # @1000K W/m/K
k_crofer=(24-23)/(800-400)*((1000-273.15)-800)+24 # @1000K W/m/K
print('porosity of the fuel channel: {:0.4g}'.format(efc))
print('porosity of the air channel: {:0.4g}'.format(efc))
print('Conductivity of Nickel at 1000K: {:0.4g} W/(m K)'.format(k_Ni))
print('Conductivity of Crofer 22 APU at 1000K: {:0.4g} W/(m K)'.format(k_crofer))
print('\nConductivity of fuel at 1000K and corresponding Ni-foam value assuming series/parallel approach with A={:0.4g}'.format(A))
for x in x_compositions: # loop each line of the compositions
    # calculate conductivity of gas mixture
    phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
    viscosity_mix=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
    conductivity_mix=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
    heat_capacity_mix=heat_capacity.dot(x)
    density_mix=1e5/8.3145/T*M.dot(x)/1000
    
    k_fuel=conductivity_mix
    
    # calculate conductivity of foam by series and parallel approach
    k_Ni_foam_series=k_Ni/((1-efc)+efc*k_Ni/k_fuel)
    k_Ni_foam_parallel=(1-efc)*k_Ni+efc*k_fuel
    
    # calculate parameter a based on A
    a=(k_Ni_foam_parallel/(A*k_Ni_foam_parallel+(1-A)*k_Ni_foam_series)-1)*k_Ni_foam_series/(k_Ni_foam_parallel-k_Ni_foam_series)
    
    # calculate k_eff based on a value
    k_eff=(a/k_Ni_foam_series+(1-a)/k_Ni_foam_parallel)**-1
    
    gas_name=':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0])
    print(('{:30s}\tk_gas={:5.4g} W/(m K) \t-->>k_series={:5.4g} W/(m K), k_parallel={:5.4g} W/(m K). k_eff={:5.4g} W/(m K) <<--a={:5.4g}').format(gas_name,conductivity_mix,k_Ni_foam_series,k_Ni_foam_parallel,k_eff,a))

print('\nConductivity of air at 1000K and corresponding Crofer 22 APU foam value assuming series/parallel approach with A={:0.4g}'.format(A))
# calculate conductivity of gas mixture (air)
x=array([0,0,0,0,0,0.7812,0.2096,0.0092,0,0,0,0]) # N2,O2,Ar
phi=1/sqrt(8)*(1+outer(M,1/M))**(-1/2)*(1+outer(viscosity,1/viscosity)**(1/2)*outer(1/M,M)**(1/4))**2
viscosity_mix=sum(viscosity*x/phi.dot(x)) # semiempirical mixing rule
conductivity_mix=sum(conductivity*x/phi.dot(x)) # semiempirical mixing rule
heat_capacity_mix=heat_capacity.dot(x)
density_mix=1e5/8.3145/T*M.dot(x)/1000

k_air=conductivity_mix

# calculate conductivity of foam by series and parallel approach
k_crofer_foam_series=k_crofer/((1-eac)+eac*k_crofer/k_air)
k_crofer_foam_parallel=(1-eac)*k_crofer+efc*k_air

# calculate parameter a based on A
a=(k_crofer_foam_parallel/(A*k_crofer_foam_parallel+(1-A)*k_crofer_foam_series)-1)*k_crofer_foam_series/(k_crofer_foam_parallel-k_crofer_foam_series)

# calculate k_eff based on fixed a value
k_eff=(a/k_crofer_foam_series+(1-a)/k_crofer_foam_parallel)**-1

gas_name=':'.join([labels[i] for i in range(len(labels)) if   x[i]!=0])+' '+':'.join(['{:d}'.format(int(round(  x[i]*100,0))) for i in range(len(labels)) if round(  x[i]*100,0)!=0])
print(('{:30s}\tk_gas={:5.4g} W/(m K) \t-->>k_series={:5.4g} W/(m K), k_parallel={:5.4g} W/(m K). k_eff={:5.4g} W/(m K) <<--a={:5.4g}').format(gas_name,conductivity_mix,k_crofer_foam_series,k_crofer_foam_parallel,k_eff,a))
```

    porosity of the fuel channel: 0.93
    porosity of the air channel: 0.93
    Conductivity of Nickel at 1000K: 72 W/(m K)
    Conductivity of Crofer 22 APU at 1000K: 23.82 W/(m K)
    
    Conductivity of fuel at 1000K and corresponding Ni-foam value assuming series/parallel approach with A=0.35
    H2:CH4:CO2:CO:H2O 49:8:8:6:29 	k_gas=0.2051 W/(m K) 	-->>k_series=0.2205 W/(m K), k_parallel=5.231 W/(m K). k_eff=1.974 W/(m K) <<--a=0.07261
    
    Conductivity of air at 1000K and corresponding Crofer 22 APU foam value assuming series/parallel approach with A=0.35
    N2:O2:Ar 78:21:1              	k_gas=0.06083 W/(m K) 	-->>k_series=0.0654 W/(m K), k_parallel=1.724 W/(m K). k_eff=0.6458 W/(m K) <<--a=0.06582


# References

[1] Bird R. B., Stweart W. E., Lightfoot E. N. (2002). Transport phenomena. 2nd ed. John Wiley & Sons. New York. S. 26, 276, 864

[2] Poling B. E., Prausnitz J. M., O'Connell J. P. (2001). The properties of gases and liquids. 5th ed. McGraw-Hill. New York. S. B.1, 10.4

[3] Burcat, A., & Ruscic, B. (2001). Third millennium ideal gas and condensed phase thermochemical database for combustion. Technion-Israel Institute of Technology.

[4] Stephan, P. u. a. (2013). VDI-Wärmeatlas. 11. bearb. und erw. Aufl. Berlin, Heidelberg: Springer-Verlag. S. 652.

[5] Perry, R. H. Green, D. W., Ed. (2008). Perry's Chemical Engineers' Handbook, Eighth Edition McGraw-Hill New York. S. 2-460.

[6] Bhattacharya, A., Calmidi, V. V., & Mahajan, R. L. (2002). Thermophysical properties of high porosity metal foams

[8] Shin, Y. C., Hashimoto, S. I., Yashiro, K., Amezawa, K., & Kawada, T. (2016). Thermal properties of perovskite-type oxides La0. 6Sr0. 4Co1-XFexO3-δ (0≤ x≤ 1.0). ECS Transactions, 72(7), 105.

[9] Pierre Villars (Chief Editor), PAULING FILE in: Inorganic Solid Phases, SpringerMaterials (online database), Springer, Heidelberg (ed.) SpringerMaterials La0.6Sr0.4Co0.2Fe0.8O3 (Sr0.6La0.4Fe0.4Co0.6O3 rt) Crystal Structure sd_1227605 (Springer-Verlag GmbH, Heidelberg, © 2016)

[10] Xu, Q., Huang, D. P., Chen, W., Zhang, F., & Wang, B. T. (2007). Structure, electrical conducting and thermal expansion properties of Ln0. 6Sr0. 4Co0. 2Fe0. 8O3 (Ln= La, Pr, Nd, Sm) perovskite-type complex oxides. Journal of alloys and compounds, 429(1-2), 34-39.

[11] M. Radovic, E. Lara/Curzio, R. Trejo, H. Wang and W. D. Porter, “Thermo/Physical Properties of Ni/YSZ as a Function of Temperature and Porosity,” Ceramic Engineering and Science Proceedings, 27, Issue 4 (2006) pp. 79/85

[13] Capdevila, X. G., Folch, J., Calleja, A., Llorens, J., Segarra, M., Espiell, F., & Morante, J. R. (2009). High-density YSZ tapes fabricated via the multi-folding lamination process. Ceramics International, 35(3), 1219-1226.

[14] Talebi, T., Sarrafi, M. H., Haji, M., Raissi, B., & Maghsoudipour, A. (2010). Investigation on microstructures of NiO–YSZ composite and Ni–YSZ cermet for SOFCs. International Journal of Hydrogen Energy, 35(17), 9440-9447.

[15] Kawashima, T., & Hishinuma, M. (1996). Thermal properties of porous Ni/YSZ particulate composites at high temperatures. Materials transactions, JIM, 37(9), 1518-1524.

[18] Yoshito, W., Matos, J., Ussui, V., Lazar, D., & Paschoal, J. (2009). Reduction kinetics of NiO–YSZ composite for application in solid oxide fuel cell. Journal of thermal analysis and calorimetry, 97(1), 303-308.

[19] N. Kostka and K. Hammeke, Ermittlung van Stoffdaten der Gase O2, N2, CO, CO2, He, H2O, H2, CH4, CH3OH und C2H5OH und der Mischungen aus diesen Gasen, JüL-Rep. 2112, Research Centre Jülich, 1987. 


```python

```

6f48516ac3707c0fcd53c3e73a6384efadd5a5ba
ipynb
Jupyter Notebook
Ni-foam.ipynb
santiago-salas-v/walas
65c1a86f9f2aaf6757d139aa6ce658ab1c750f0d
MIT
Ni-foam.ipynb
santiago-salas-v/walas
65c1a86f9f2aaf6757d139aa6ce658ab1c750f0d
MIT
Ni-foam.ipynb
santiago-salas-v/walas
65c1a86f9f2aaf6757d139aa6ce658ab1c750f0d
MIT
Qwen/Qwen-72B
1. YES
2. YES