# Demo entry 6636439

11111

Submitted by anonymous on Aug 24, 2017 at 14:32
Language: Python. Code size: 5.7 kB.

```#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Sun May 21 14:01:04 2017
@author: LINWEI
"""

import numpy as np
from scipy.optimize import brentq
from scipy.optimize import newton_krylov
import matplotlib.pyplot as plt

# Global Constant
epsilon0 = 8.85E-12  # dielectric permittiviy of vacuum
epsilonr = 81.  # relative dielectric permittivity of water
k = 1.38E-23  # Boltzmann constant
T = 298.15  # temperature
eta = 0.8949E-3  # viscosity of water
ec = 1.602E-19  # elementary charge
NA = 6.022E23  # Avogadro constant

def phi_func(a, b):
"""
Parameters
----------
a ：hydrodynamic radius of the particle
b : radius of the Wigner-Seitz spherical cell

Return
------
phi : volume fraction of particles in the suspension
"""
return (a/b)**3

def b_func(a, phi):
"""
Parameters
----------
a ：hydrodynamic radius of the particle
phi : volume fraction of particles in the suspension

Return
------
b : radius of the Wigner-Seitz spherical cell
"""
return a/phi**(1.0/3.0)

def theory(a, Z, phi):
sigma = -Z*ec/(4.*np.pi*a**2)
b = b_func(a, phi)  # radius of the Wigner-Seitz spherical cells
n_P = phi/((4./3.)*np.pi*a**3)  # number concentration of particles
item = (Z*n_P)**2 + 4E-8*NA**2
item = np.sqrt(item)
n_bar_H = (Z*n_P + item)*0.5    # number concentration of hydrogen ions
n_bar_OH = (-Z*n_P + item)*0.5    # number concentratioon of hydroxide ions
N_bar_H = n_bar_H/n_bar_OH

m = 500  # mesh number
RLeft, RRight = 1., b/a  # endpoints
yPrimeLeft, yPrimeRight = -(sigma*ec*a)/(epsilon0*epsilonr*k*T), 0.
h = (RRight-RLeft)/m
R = np.arange(RLeft, RRight+h, h)

def trapezoid(x, y):
"""
Composite Trapezoidal rule for integration, Ingegrate_x[0]^x[-1] (y(x)dx)
"""
return 0.5*h*(2.*y.sum()-y[0]-y[-1])

def Psi_func(N_Psi0_H):
def F_PB(x, y, yPrime):
item1 = -2.0*yPrime/x
kappa2 = (ec**2*n_bar_OH)/(epsilon0*epsilonr*k*T)
item2 = -kappa2*a**2
item3 = N_Psi0_H*np.exp(-y) - np.exp(y)
return item1+item2*item3

def residual_PB(y):
r = np.zeros_like(y)
r[0] = 2.*(y[1]-y[0]) - h**2*F_PB(R[0], y[0], yPrimeLeft) - 2.*h*yPrimeLeft
r[-1] = 2.*(y[-2]-y[-1])-h**2*F_PB(R[-1], y[-1], yPrimeRight) + 2.*h*yPrimeRight
r[1:-1] = (y[2:]-2.*y[1:-1]+y[:-2]) - h**2*F_PB(R[1:-1], y[1:-1], (y[2:]-y[:-2])/(2.*h))
return r

def start_PB(x):
y = -5. + 5./(RRight-RLeft)*(x-RLeft)
return y

y = newton_krylov(residual_PB, start_PB(R))
return y

def residual_N(N_Psi0_H):
def f(x, y):
return (N_Psi0_H*np.exp(-y)-np.exp(y))*x**2

r = 4.*np.pi*a**3*n_bar_OH*trapezoid(R, f(R, Psi_func(N_Psi0_H))) - Z
return r

def findRange():
N_Psi0_H = N_bar_H*0.3
zla = residual_N(N_Psi0_H)
while np.sign(zla) == 1:
N_Psi0_H -= N_bar_H*0.1
zla = residual_N(N_Psi0_H)
if np.sign(zla) == -1:
break
N_Psi0_H_1 = N_Psi0_H
N_Psi0_H_2 = N_Psi0_H + N_bar_H*0.1
return N_Psi0_H_1, N_Psi0_H_2

N_Psi0_H_1, N_Psi0_H_2 = findRange()
N_Psi0_H = brentq(residual_N, N_Psi0_H_1, N_Psi0_H_2)
print "N_Psi0_H = ", N_Psi0_H
Psi = Psi_func(N_Psi0_H)

# Rho[R]
Rho = N_Psi0_H*np.exp(-Psi)-np.exp(Psi)
# Rho'[R]
RhoPrime = np.zeros_like(Rho)
RhoPrime[0] = (-3.*Rho[0]+4.*Rho[1]-Rho[2])/(2.*h)
RhoPrime[-1] = (-3.*Rho[-1]+4.*Rho[-2]-Rho[-3])/(2.*h)
RhoPrime[1:-1] = (Rho[2:]-Rho[:-2])/(2.*h)
# Gamma[R]
Gamma = (1./(1.-phi))*(1.+0.5/R**3)*RhoPrime
integral_1 = np.zeros_like(R)
for i in range(0, R.size):
integral_1[i] = trapezoid(R[i:], Gamma[i:])
integral_2 = np.zeros_like(R)
GammaR3 = Gamma*R**3
for i in range(0, R.size):
integral_2[i] = trapezoid(R[i:], GammaR3[i:])
integral = -trapezoid(R, (R-R**4*phi)*integral_1)+trapezoid(R, (1./R**2-R*phi)*integral_2)
mu1 = -(2./9.)*(ec*n_bar_OH*a**2/eta)*integral
rhob = Rho[-1]*ec*n_bar_OH
Yb = (1./(1.-phi))*(b + 0.5*a**3/b**2)
mu2 = -(2./9.)*((rhob*Yb)/(eta*b**2))*(a**5/(5.*b**2)-(a**2*b)+1.8*b**3-b**4/a)
mu = mu1 + mu2  # electric mobility
return mu

if __name__ == "__main__":
a = 42.5E-9  # hydrodynamic radius
Z = 895.
phi = 0.001
mu = theory(a, Z, phi)
print mu
phi2 = np.arange(0.001, 0.01, 0.001)
phi3 = np.arange(0.01, 0.06, 0.01)
phi_range = np.concatenate((phi2, phi3))
mu_range = np.zeros_like(phi_range)

# --------- #
# my theory #
# --------- #
for i, phi in enumerate(phi_range):
print
print 'phi = ', phi
mu_range[i] = theory(a, Z, phi)
print 'mu = ', mu_range[i]

plt.semilogx(phi_range, mu_range, label="theory")

# -------------------- #
# experimental results #
# -------------------- #
phi_experiment = np.array([0.002, 0.003, 0.004, 0.005, 0.006,
0.008, 0.009, 0.02, 0.03, 0.04,
0.05])
mu_experiment = np.array([10.13, 9.41, 8.97, 9.25, 7.26,
7.83, 6.93, 7.63, 6.82, 6.84,
5.87])
mu_experiment = 1E-8*mu_experiment
plt.semilogx(phi_experiment, mu_experiment, 'o')
plt.legend()

plt.show()
```

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.