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.

Delete this entry (admin only).