contact-angle.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy import init_printing\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "from scipy.optimize import newton\n",
    "\n",
    "init_printing(use_unicode=False, wrap_line=False, no_global=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calculate_helmholtz_excess(\n",
    "    rho_value,     \n",
    "    reduced_temperature, \n",
    "    rho_inf, \n",
    "    eos\n",
    "    ):\n",
    "    # Defining symbols\n",
    "    R, T, v, b, a, alpha, rho = symbols(\"R T v b a alpha rho\")\n",
    "    M, k = symbols(\"M xi\")\n",
    "    P_r, rho_r, T_r = symbols(\"P_r rho_r T_r\")\n",
    "    P_c, rho_c, T_c = symbols(\"P_c rho_c T_c\")\n",
    "    omega_a, omega_b = symbols(\"Omega_a Omega_b\")\n",
    "    Z_c = symbols(\"Z_c\")\n",
    "    rho_oo = symbols(\"rho_oo\")\n",
    "    \n",
    "    P1 = eos.pressure.subs(v, M/rho)\n",
    "    P2 = P1.subs(rho, rho_c * rho_r)\n",
    "    P3 = P2.subs(T, T_c * T_r)\n",
    "    P4 = P3/P_c\n",
    "    P5 = P4.subs(a, omega_a*(R**2)*(T_c**2)/P_c)\n",
    "    P6 = P5.subs(b, omega_b*R*T_c/P_c)\n",
    "    P7 = expand(P6)\n",
    "    P8 = P7.subs(T_c, P_c * M / (Z_c * rho_c * R)) \n",
    "    \n",
    "    P_r = simplify(P8)    \n",
    "    P_oo = P_r.subs(rho_r, rho_oo)    \n",
    "    integrand = ((P_r - P_oo)/(rho_r ** 2))\n",
    "\n",
    "    integrand = integrand.subs(omega_a, eos.omega_a)\n",
    "    integrand = integrand.subs(omega_b, eos.omega_b)\n",
    "    integrand = integrand.subs(Z_c, eos.critical_z_factor)\n",
    "    integrand = integrand.subs(T_r, reduced_temperature)\n",
    "    integrand = integrand.subs(rho_oo, rho_inf)\n",
    "    \n",
    "    func = lambdify(rho_r, integrand)\n",
    "    y, err = integrate.quad(func, rho_inf, rho_value)\n",
    "    return rho_value * y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Defining symbols\n",
    "P, R, T, v, b, a, alpha, rho = symbols(\"P R T v b a alpha rho\")\n",
    "def find_correct_root_of_cubic_eos(p0, p1, p2, p3, fluid_type):\n",
    "    coef_a = (3.0 * p2 - (p1 ** 2)) / 3.0        \n",
    "    coef_b = (2.0 * (p1 ** 3) - 9.0 * p1 * p2 + 27.0 * p3) / 27.0        \n",
    "    delta = 0.25 * (coef_b ** 2) + (coef_a ** 3) / 27.0     \n",
    "\n",
    "    if delta > 0.0:\n",
    "        # 1 real root, 2 imaginary                 \n",
    "        const_A =  cbrt(-0.5 * coef_b + sqrt(delta)) \n",
    "        const_B =  cbrt(-0.5 * coef_b - sqrt(delta))\n",
    "\n",
    "        correct_root = const_A + const_B - p1 / 3.0 \n",
    "    else:\n",
    "        # 3 real roots\n",
    "        phi = acos(-0.5 * coef_b / sqrt(-(coef_a ** 3) / 27.0))\n",
    "        root_1 = 2.0 * sqrt(-coef_a / 3.0) * cos(phi / 3.0) - p1 / 3.0\n",
    "        root_2 = 2.0 * sqrt(-coef_a / 3.0) * cos(phi / 3.0 + 2.0 * np.pi / 3.0) - p1 / 3.0\n",
    "        root_3 = 2.0 * sqrt(-coef_a / 3.0) * cos(phi / 3.0 + 4.0 * np.pi / 3.0) - p1 / 3.0\n",
    "\n",
    "        smallest_root = min(min(root_1,root_2), root_3)\n",
    "        largest_root = max(max(root_1,root_2), root_3)\n",
    "\n",
    "        if fluid_type is 'liquid':        \n",
    "            correct_root = smallest_root\n",
    "        else:\n",
    "            assert fluid_type is 'vapor', 'Wrong fluid type! ' + fluid_type\n",
    "            correct_root = largest_root\n",
    "    \n",
    "    assert correct_root > 0.0, fluid_type + ' Z-factor < 0.0! Delta is %f, %f' % (delta, correct_root)\n",
    "    \n",
    "    return correct_root\n",
    "\n",
    "GAS_CONST = 8.314\n",
    "\n",
    "class EquationOfState(object):\n",
    "    def __init__(\n",
    "        self,\n",
    "        critical_pressure, \n",
    "        critical_temperature,\n",
    "        acentric_factor,\n",
    "        molar_mass\n",
    "    ):       \n",
    "        self.critical_pressure = critical_pressure\n",
    "        self.critical_temperature = critical_temperature\n",
    "        self.acentric_factor = acentric_factor\n",
    "        self.molar_mass = molar_mass\n",
    "        \n",
    "        \n",
    "    def calculate_pressure_and_densities(self, temperature):\n",
    "        \n",
    "        def func(x):            \n",
    "            pressure = x\n",
    "            Z_V = eos.calculate_eos_roots(pressure, temperature, 'vapor')\n",
    "            Z_L = eos.calculate_eos_roots(pressure, temperature, 'liquid')\n",
    "        \n",
    "            f_V = eos.calculate_fugacity(pressure, temperature, Z_V)\n",
    "            f_L = eos.calculate_fugacity(pressure, temperature, Z_L)\n",
    "            \n",
    "            return f_V - f_L\n",
    "        \n",
    "        pressure = newton(func=func, x0=1.0e5)\n",
    "        \n",
    "        Z_V = eos.calculate_eos_roots(pressure, temperature, 'vapor')\n",
    "        Z_L = eos.calculate_eos_roots(pressure, temperature, 'liquid')\n",
    "        \n",
    "        rho_V = pressure * self.molar_mass / (Z_V * GAS_CONST * temperature)\n",
    "        rho_L = pressure * self.molar_mass / (Z_L * GAS_CONST * temperature)\n",
    "        \n",
    "        return pressure, rho_V, rho_L\n",
    "    \n",
    "    \n",
    "class PengRobinsonEos(EquationOfState):\n",
    "    def __init__(\n",
    "        self,\n",
    "        critical_pressure, \n",
    "        critical_temperature,\n",
    "        acentric_factor,\n",
    "        molar_mass\n",
    "    ):\n",
    "        self.omega_a = 0.45724\n",
    "        self.omega_b = 0.07780\n",
    "        self.critical_z_factor = 0.3074\n",
    "        self.pressure = R * T / (v - b) - a*alpha / (v*(v + b) + b*(v - b))\n",
    "        \n",
    "        EquationOfState.__init__(self, critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "\n",
    "    \n",
    "    def alpha_function(self, reduced_temperature):\n",
    "        omega = self.acentric_factor\n",
    "        \n",
    "        if omega < 0.49:\n",
    "            m = 0.37464 + 1.54226 * omega - 0.26992 * (omega ** 2)\n",
    "        else:\n",
    "            m = 0.379642 + 1.48503 * omega - 0.164423 * (omega ** 2) + 0.016667 * (omega ** 3)\n",
    "        \n",
    "        return (1.0 + m * (1.0 - sqrt(reduced_temperature))) ** 2\n",
    "    \n",
    "    \n",
    "    def calculate_fugacity(\n",
    "        self,\n",
    "        pressure,\n",
    "        temperature,\n",
    "        z_factor\n",
    "    ):\n",
    "        P = pressure\n",
    "        T = temperature\n",
    "        Z = z_factor\n",
    "        \n",
    "        Pr = P / self.critical_pressure\n",
    "        Tr = T / self.critical_temperature\n",
    "        \n",
    "        aa = self.alpha_function(Tr)\n",
    "        \n",
    "        A = aa * self.omega_a * Pr / (Tr ** 2)\n",
    "        B = self.omega_b * Pr / Tr\n",
    "        \n",
    "        SQRT_2 = sqrt(2.0)        \n",
    "        ln_f = (Z - 1.0) - log( Z - B ) - (A / (2.0 * SQRT_2 * B)) \\\n",
    "             * log( (Z + (1.0 + SQRT_2)*B)/(Z + (1.0 - SQRT_2)*B) )\n",
    "        return P * exp(ln_f)\n",
    "        \n",
    "        \n",
    "    def calculate_eos_roots(self, pressure, temperature, fluid_type):\n",
    "        P = pressure\n",
    "        T = temperature\n",
    "\n",
    "        Pr = P / self.critical_pressure\n",
    "        Tr = T / self.critical_temperature\n",
    "        \n",
    "        aa = self.alpha_function(Tr)\n",
    "        \n",
    "        A = aa * self.omega_a * Pr / (Tr ** 2)\n",
    "        B = self.omega_b * Pr / Tr      \n",
    "        \n",
    "        p0 = 1.0\n",
    "        p1 = - (1.0 - B)\n",
    "        p2 = A - 3.0 * (B ** 2) - 2.0 * B\n",
    "        p3 = -(A * B - B ** 2 - B ** 3)     \n",
    "\n",
    "        return find_correct_root_of_cubic_eos(p0, p1, p2, p3, fluid_type)       \n",
    "        \n",
    "        \n",
    "class VanDerWaalsEos(EquationOfState):\n",
    "    def __init__(\n",
    "        self,\n",
    "        critical_pressure, \n",
    "        critical_temperature,\n",
    "        acentric_factor,\n",
    "        molar_mass\n",
    "    ):\n",
    "        self.omega_a = 27.0/64.0\n",
    "        self.omega_b = 1.0/8.0\n",
    "        self.critical_z_factor = 3.0/8.0\n",
    "        self.pressure = R * T / (v - b) - a / (v**2)\n",
    "       \n",
    "        EquationOfState.__init__(self, critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "        \n",
    "    def calculate_fugacity(\n",
    "        self,\n",
    "        pressure,\n",
    "        temperature,\n",
    "        z_factor\n",
    "    ):\n",
    "        P = pressure\n",
    "        T = temperature\n",
    "        Z = z_factor\n",
    "        Pr = P / self.critical_pressure\n",
    "        Tr = T / self.critical_temperature\n",
    "        \n",
    "        A = self.omega_a * Pr / (Tr ** 2)\n",
    "        B = self.omega_b * Pr / Tr\n",
    "        \n",
    "        ln_f = B / (Z - B) - log(Z - B) - 2.0 * A / Z\n",
    "        return P * exp(ln_f)\n",
    "    \n",
    "    \n",
    "    def calculate_eos_roots(self, pressure, temperature, fluid_type):\n",
    "        P = pressure\n",
    "        T = temperature\n",
    "\n",
    "        Pr = P / self.critical_pressure\n",
    "        Tr = T / self.critical_temperature\n",
    "        \n",
    "        A = self.omega_a * Pr / (Tr ** 2)\n",
    "        B = self.omega_b * Pr / Tr\n",
    "        \n",
    "        p0 = 1.0\n",
    "        p1 = - (B + 1.0)\n",
    "        p2 = A\n",
    "        p3 = - A * B   \n",
    "        \n",
    "        return find_correct_root_of_cubic_eos(p0, p1, p2, p3, fluid_type) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "temperature = 303.15\n",
    "pressure = 1.0 * 1.0e5\n",
    "\n",
    "# n-Butane properties\n",
    "critical_pressure = 6894.75729 * 550.7\n",
    "critical_temperature = (5.0 / 9.0) * 765.3 # [K]\n",
    "acentric_factor = 0.192800 # [-]\n",
    "molar_mass = 58.12 * 0.001 # [kg/mol]\n",
    "\n",
    "eos = VanDerWaalsEos(critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "\n",
    "pressure, rho_V, rho_L = eos.calculate_pressure_and_densities(temperature)\n",
    "Z_V = eos.calculate_eos_roots(pressure, temperature, 'vapor')\n",
    "Z_L = eos.calculate_eos_roots(pressure, temperature, 'liquid')\n",
    "f_V = eos.calculate_fugacity(pressure, temperature, Z_V)\n",
    "f_L = eos.calculate_fugacity(pressure, temperature, Z_L)\n",
    "\n",
    "assert abs(f_L - f_V) < 1.0e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step by step\n",
    "\n",
    "1. Given the temperature and an Equation of State\n",
    "2. Calculate the pressure and the densities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Calculate Helmholtz Excess Energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# n-Butane properties\n",
    "critical_pressure = 6894.75729 * 550.7\n",
    "critical_temperature = (5.0 / 9.0) * 765.3 # [K]\n",
    "acentric_factor = 0.192800 # [-]\n",
    "molar_mass = 58.12 * 0.001 # [kg/mol]\n",
    "\n",
    "reduced_temperatures = np.array([0.8, 0.82, 0.87, 0.9])\n",
    "\n",
    "for reduced_temperature in reduced_temperatures:\n",
    "    temperature = reduced_temperature * critical_temperature    \n",
    "    \n",
    "    eos = PengRobinsonEos(critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "    pressure, rho_V, rho_L = eos.calculate_pressure_and_densities(temperature)\n",
    "\n",
    "    rho_c = eos.critical_pressure * eos.molar_mass / (eos.critical_z_factor * GAS_CONST * eos.critical_temperature)\n",
    "    rho_inf = rho_V / rho_c \n",
    "    \n",
    "    eos.pressure = eos.pressure.subs(alpha, eos.alpha_function(reduced_temperature))    \n",
    "    \n",
    "    reduced_densities = np.linspace(0.1, 2.6, num=40)\n",
    "    \n",
    "    psi_vec = []\n",
    "    reduced_densities_to_plot = []\n",
    "    \n",
    "    for rho in reduced_densities:\n",
    "        psi = calculate_helmholtz_excess(\n",
    "            rho_value=rho, \n",
    "            reduced_temperature=reduced_temperature, \n",
    "            rho_inf=rho_inf,\n",
    "            eos=eos\n",
    "            )\n",
    "\n",
    "        reduced_densities_to_plot.append(rho)\n",
    "        psi_vec.append(psi)\n",
    "        \n",
    "    psi_vec = np.array(psi_vec)\n",
    "    plt.plot(reduced_densities_to_plot, psi_vec, marker='o', \n",
    "             label='Tr = %.2f' %reduced_temperature)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.grid(True)\n",
    "ax.set_xlabel(xlabel='Reduced Density [-]')\n",
    "ax.set_ylabel(ylabel=r'$\\Psi^*_{exc}$')\n",
    "ax.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def integrate_helholtz(eos, reduced_temperature, rho_0_r, rho_1_r, rho_inf):\n",
    "    def func(x):\n",
    "        psi = calculate_helmholtz_excess(\n",
    "                rho_value=x, \n",
    "                reduced_temperature=reduced_temperature, \n",
    "                rho_inf=rho_inf,\n",
    "                eos=eos\n",
    "                )\n",
    "        return psi\n",
    "        \n",
    "    y, err = integrate.quad(func, rho_0_r, rho_1_r)\n",
    "    return sqrt(2.0 * y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surface parameter $\\kappa$ in terms of the surface tension\n",
    "\n",
    "$$\n",
    "\\gamma_{lv} = \\rho_c \\sqrt{\\kappa P_c} \\int^{\\rho^*_{l}}_{\\rho^*_{v}} \\sqrt{2\\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "We would like to correlate the surface parameter with experimental values of surface tension:\n",
    "$$\n",
    "\\kappa = \\left(\\frac{\\gamma_{lv,exp}}{\\rho_c \\sqrt{P_c} \\int^{\\rho^*_{l}}_{\\rho^*_{v}} \\sqrt{2\\Psi^*_{exc}} d\\rho^*} \\right)^2\n",
    "$$\n",
    "\n",
    "where $\\gamma_{lv,exp}$ is the experimental value of the surface tension.\n",
    "\n",
    "Since I was not able to find experimental values of surface tension of n-butane for more than one temperature point, the component used will be water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Experimental surface tension n-butane\n",
    "# T [K]   gamma [N/m]\n",
    "# 230       0.02064\n",
    "# 288       0.012\n",
    "#reduced_temperature_points = np.array([230.0, 288]) / critical_temperature\n",
    "#surface_tension_points = np.array([0.02064, 0.012])\n",
    "\n",
    "# Experimental surface tension Water - H2O\n",
    "reduced_temperature_points = np.array([0.422, 0.4297, 0.4374, 0.4529, 0.4683, \n",
    "                                       0.4838, 0.4992, 0.5147, 0.5301, 0.5456, \n",
    "                                       0.5610, 0.5765, 0.6537, 0.7310, 0.8082,\n",
    "                                      0.8854, 0.9627, 0.9999])\n",
    "surface_tension_experimental_points = np.array([0.0756, 0.0749, 0.0742, 0.0728, 0.0712,\n",
    "                                   0.0696, 0.0679, 0.0662, 0.0644, 0.0626, \n",
    "                                   0.0608, 0.0589, 0.0482, 0.0376, 0.0264, \n",
    "                                   0.0147, 0.0037, 0.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# H2O properties\n",
    "critical_pressure = 221.2 * 1.0e5\n",
    "critical_temperature = 647.3 # [K]\n",
    "acentric_factor = 0.344 # [-]\n",
    "molar_mass = 18.0 * 0.001 # [kg/mol]\n",
    "\n",
    "reduced_temperature_values = reduced_temperature_points[3:-2] # Two-phase values only\n",
    "surface_tension_points = surface_tension_experimental_points[3:-2] # Two-phase values only\n",
    "\n",
    "#reduced_temperature_values = reduced_temperature_points[3:4] # Two-phase values only\n",
    "#surface_tension_points = surface_tension_experimental_points[3:4] # Two-phase values only\n",
    "\n",
    "kappa_values = []\n",
    "\n",
    "for reduced_temperature, surface_tension in zip(reduced_temperature_values, surface_tension_points):\n",
    "    temperature = reduced_temperature * critical_temperature    \n",
    "    \n",
    "    eos = PengRobinsonEos(critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "    eos.pressure = eos.pressure.subs(alpha, eos.alpha_function(reduced_temperature))  \n",
    "    #eos = VanDerWaalsEos(critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "    \n",
    "    pressure, rho_V, rho_L = eos.calculate_pressure_and_densities(temperature)\n",
    "\n",
    "    rho_c = eos.critical_pressure * eos.molar_mass / (eos.critical_z_factor * GAS_CONST * eos.critical_temperature)\n",
    "    #rho_c = 197.29\n",
    "    rho_V_r = rho_V / rho_c\n",
    "    rho_L_r = rho_L / rho_c\n",
    "    \n",
    "    integral = integrate_helholtz(eos, reduced_temperature, rho_V_r, rho_L_r, rho_inf=rho_V_r)\n",
    "\n",
    "    kappa = ( surface_tension / (rho_c * sqrt(eos.critical_pressure) * integral) ) ** 2\n",
    "    kappa_values.append(kappa)\n",
    "    \n",
    "kappa_values = np.array(kappa_values)    \n",
    "plt.semilogy(surface_tension_points, kappa_values, marker='o', label='Peng-Robinson')\n",
    "\n",
    "#plt.plot(gamma_calc, kappa_values, marker='o')\n",
    "ax = plt.gca()\n",
    "ax.grid(True, which=u'both')\n",
    "ax.set_xlabel(xlabel='Surface Tension [N/m]')\n",
    "ax.set_ylabel(ylabel=r'Surface parameter $\\kappa$ ($N / \\left( kg/m^3 \\right)^2$)')\n",
    "ax.legend(loc='upper left')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.semilogy(surface_tension_points, kappa_values, marker='o', label='Peng-Robinson')\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.grid(True, which=u'both')\n",
    "ax.set_xlabel(xlabel='Surface Tension [N/m]')\n",
    "ax.set_ylabel(ylabel=r'Surface parameter $\\kappa$ ($N / \\left( kg/m^3 \\right)^2$)')\n",
    "ax.legend(loc='upper left')\n",
    "ax.axis([0.01, 0.08, 1.0e-17, 1.0e-14])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Expressions used\n",
    "\n",
    "Calculate the above curve for various temperatures\n",
    "\n",
    "Calculate the equilibrium contact angle as a function of $\\phi^*$\n",
    "\n",
    "$$\n",
    "\\cos \\theta = \\frac{\\gamma_{sv} - \\gamma_{sl}}{\\gamma_{lv}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\phi^* = \\frac{\\phi}{\\sqrt{\\kappa P_c}} = \\pm \\sqrt{2 \\Psi^*_{exc}(\\rho_s, T)}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\gamma_{vs} = \\Phi(\\rho_s) + \\rho_c \\sqrt{\\kappa P_c} \\int^{\\rho^*_{sv}}_{\\rho^*_{v}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\gamma_{ls} = \\Phi(\\rho_s) + \\rho_c \\sqrt{\\kappa P_c} \\int^{\\rho^*_{sl}}_{\\rho^*_{l}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Phi(\\rho_s) = - \\phi \\rho^*_{s} \\rho_c\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\phi = \\phi^* \\sqrt{\\kappa P_c}\n",
    "$$\n",
    "\n",
    "Obs: Check out fig 3.9 to find the roots of eq $\\phi^* = \\frac{\\phi}{\\sqrt{\\kappa P_c}} = \\pm \\sqrt{2 \\Psi^*_{exc}(\\rho_s, T)}$\n",
    "\n",
    "We also can get $\\kappa$ from \n",
    "\n",
    "$$\n",
    "\\gamma_{lv} = \\rho_c \\sqrt{2 \\kappa P_c} \\int^{\\rho^*_{l}}_{\\rho^*_{v}} \\sqrt{\\Psi^*_{exc}} d\\rho^*\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "First, let's plot the curve $\\sqrt{2 \\Psi^*_{exc}(\\rho, T)}$ and intersect it with $\\phi^* = 0.2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Plot of equation to be solved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sqrt_psi_vec = {}\n",
    "reduced_densities_to_plot = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# n-Butane properties\n",
    "critical_pressure = 6894.75729 * 550.7\n",
    "critical_temperature = (5.0 / 9.0) * 765.3 # [K]\n",
    "acentric_factor = 0.192800 # [-]\n",
    "molar_mass = 58.12 * 0.001 # [kg/mol]\n",
    "\n",
    "reduced_temperature = 0.8\n",
    "temperature = reduced_temperature * critical_temperature\n",
    "\n",
    "eos = PengRobinsonEos(critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "eos.pressure = eos.pressure.subs(alpha, eos.alpha_function(reduced_temperature))  \n",
    "pressure, rho_V, rho_L = eos.calculate_pressure_and_densities(temperature)\n",
    "rho_c = eos.critical_pressure * eos.molar_mass / (eos.critical_z_factor * GAS_CONST * eos.critical_temperature)\n",
    "\n",
    "rho_V_r = rho_V / rho_c\n",
    "rho_L_r = rho_L / rho_c\n",
    "\n",
    "reduced_densities = np.linspace(0.00001, 3.5, num=150)\n",
    "    \n",
    "sqrt_psi_vec[reduced_temperature] = []\n",
    "reduced_densities_to_plot[reduced_temperature] = []\n",
    "\n",
    "for rho in reduced_densities:\n",
    "    psi = calculate_helmholtz_excess(\n",
    "        rho_value=rho, \n",
    "        reduced_temperature=reduced_temperature, \n",
    "        rho_inf=rho_V_r,\n",
    "        eos=eos\n",
    "        )\n",
    "\n",
    "    reduced_densities_to_plot[reduced_temperature].append(rho)\n",
    "    sqrt_psi_vec[reduced_temperature].append(sqrt(2.0*psi))\n",
    "\n",
    "sqrt_psi_vec[reduced_temperature] = np.array(sqrt_psi_vec[reduced_temperature])\n",
    "\n",
    "plt.plot(reduced_densities_to_plot[reduced_temperature], \n",
    "         sqrt_psi_vec[reduced_temperature], marker='o', label='Reduced temperature = %.2f' %reduced_temperature)\n",
    "\n",
    "surf_energy = 0.3\n",
    "phi_r = surf_energy * np.ones(len(reduced_densities_to_plot[reduced_temperature]))\n",
    "plt.plot(reduced_densities_to_plot[reduced_temperature], phi_r, 'k--', label='Specific Surface Energy = %.2f' %surf_energy)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.grid(True)\n",
    "ax.set_xlabel(xlabel='Reduced density [-]')\n",
    "ax.set_ylabel(ylabel=r'$\\sqrt{2\\Psi^*_{exc}}$')\n",
    "ax.legend(loc='upper left')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to find the roots of the equation $\\phi^* - \\sqrt{2 \\Psi^*_{exc}(\\rho, T)} = 0$ and find $\\rho_s$. So let's define a residual equation as $R(\\rho) = \\phi^* - \\sqrt{2 \\Psi^*_{exc}(\\rho, T)}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def calculate_surface_densities(phi_r, reduced_temperature, rho_V_r, rho_L_r, rho_inf):\n",
    "    def func(x):\n",
    "        psi = calculate_helmholtz_excess(\n",
    "            rho_value=x, \n",
    "            reduced_temperature=reduced_temperature, \n",
    "            rho_inf=rho_inf,\n",
    "            eos=eos\n",
    "            )\n",
    "        \n",
    "        if phi_r > 0.0:\n",
    "            return phi_r - sqrt(2.0 * psi)\n",
    "        else:\n",
    "            return phi_r + sqrt(2.0 * psi)\n",
    "\n",
    "    if phi_r > 0.0:\n",
    "        rho_SV_r = newton(func=func, x0=rho_V_r + 0.1*rho_V_r)\n",
    "        rho_SL_r = newton(func=func, x0=rho_L_r + 0.1*rho_L_r)\n",
    "    else:       \n",
    "        from scipy.optimize import brentq\n",
    "        rho_SV_r = brentq(f=func, a=0.000001, b=0.9999*rho_V_r)\n",
    "        rho_SL_r = brentq(f=func, a=0.5*(rho_L_r + rho_V_r), b=0.9999*rho_L_r)\n",
    "    \n",
    "    return rho_SV_r, rho_SL_r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need\n",
    "$$\n",
    "\\gamma_{vs} = \\Phi(\\rho_s) + \\rho_c \\sqrt{\\kappa P_c} \\int^{\\rho^*_{sv}}_{\\rho^*_{v}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\gamma_{ls} = \\Phi(\\rho_s) + \\rho_c \\sqrt{\\kappa P_c} \\int^{\\rho^*_{sl}}_{\\rho^*_{l}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Phi(\\rho_s) = - \\phi \\rho^*_{s} \\rho_c = - \\rho^*_{s} \\rho_c \\phi^* \\sqrt{\\kappa P_c}\n",
    "$$\n",
    "$$\n",
    "\\gamma_{lv} = \\rho_c \\sqrt{2 \\kappa P_c} \\int^{\\rho^*_{l}}_{\\rho^*_{v}} \\sqrt{\\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "$$\n",
    "\\cos \\theta = \\frac{\\gamma_{sv} - \\gamma_{sl}}{\\gamma_{lv}}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\gamma_{vs}^* = \\frac{\\gamma_{vs}}{\\rho_c\\sqrt{\\kappa P_c}} = - \\rho^*_{s} \\phi^* + \\int^{\\rho^*_{sv}}_{\\rho^*_{v}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\gamma_{ls}^* = \\frac{\\gamma_{ls}}{\\rho_c\\sqrt{\\kappa P_c}} =- \\rho^*_{s} \\phi^* + \\int^{\\rho^*_{sl}}_{\\rho^*_{l}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\gamma_{lv}^* = \\frac{\\gamma_{lv}}{\\rho_c\\sqrt{\\kappa P_c}} = \\int^{\\rho^*_{l}}_{\\rho^*_{v}} \\sqrt{2 \\Psi^*_{exc}} d\\rho^*\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\cos \\theta = \\frac{\\gamma_{sv} - \\gamma_{sl}}{\\gamma_{lv}} = \\frac{\\gamma_{vs}^* - \\gamma_{ls}^*}{\\gamma_{lv}^*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calculate_contact_angle(phi_r, gamma_VL_r, eos, reduced_temperature, rho_V_r, rho_L_r):\n",
    "    if abs(phi_r) < 1.0e-8:\n",
    "        contact_angle = 90.0\n",
    "    else:\n",
    "        rho_SV_r, rho_SL_r = calculate_surface_densities(\n",
    "            phi_r=phi_r, \n",
    "            reduced_temperature=reduced_temperature, \n",
    "            rho_V_r=rho_V_r, \n",
    "            rho_L_r=rho_L_r, \n",
    "            rho_inf=rho_V_r\n",
    "        )\n",
    "\n",
    "        if phi_r > 0.0:\n",
    "            integral_SV = integrate_helholtz(eos, reduced_temperature, rho_V_r, rho_SV_r, rho_inf=rho_V_r)\n",
    "            integral_SL = integrate_helholtz(eos, reduced_temperature, rho_L_r, rho_SL_r, rho_inf=rho_V_r)\n",
    "        else:\n",
    "            integral_SV = integrate_helholtz(eos, reduced_temperature, rho_SV_r, rho_V_r, rho_inf=rho_V_r)\n",
    "            integral_SL = integrate_helholtz(eos, reduced_temperature, rho_SL_r, rho_L_r, rho_inf=rho_V_r)    \n",
    "\n",
    "        gamma_SV_r = integral_SV - phi_r * rho_SV_r\n",
    "        gamma_SL_r = integral_SL - phi_r * rho_SL_r\n",
    "\n",
    "        contact_angle = acos((gamma_SV_r - gamma_SL_r) / gamma_VL_r) * 180.0 / np.pi\n",
    "    \n",
    "    return contact_angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Contact angle as a function of the interfacial tension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "contact_angles = {}\n",
    "phi_values = {}\n",
    "phi_limits = {}\n",
    "# Phi limits (so functions do not diverge)\n",
    "# Tr = 0.9 -> phi_r = 0.35\n",
    "# Tr = 0.87 -> phi_r = 0.43\n",
    "# Tr = 0.82 -> phi_r = 0.565\n",
    "# Tr = 0.8 -> phi_r = 0.615\n",
    "phi_limits[0.9] = 0.35\n",
    "phi_limits[0.87] = 0.43\n",
    "phi_limits[0.82] = 0.565\n",
    "phi_limits[0.8] = 0.615"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# n-Butane properties\n",
    "critical_pressure = 6894.75729 * 550.7\n",
    "critical_temperature = (5.0 / 9.0) * 765.3 # [K]\n",
    "acentric_factor = 0.192800 # [-]\n",
    "molar_mass = 58.12 * 0.001 # [kg/mol]\n",
    "\n",
    "\n",
    "reduced_temperature_values = [0.8, 0.82, 0.87, 0.9]\n",
    "for reduced_temperature in reduced_temperature_values:\n",
    "    temperature = reduced_temperature * critical_temperature\n",
    "    eos = PengRobinsonEos(critical_pressure, critical_temperature, acentric_factor, molar_mass)\n",
    "    rho_c = eos.critical_pressure * eos.molar_mass / (eos.critical_z_factor * GAS_CONST * eos.critical_temperature)\n",
    "\n",
    "    pressure, rho_V, rho_L = eos.calculate_pressure_and_densities(temperature)\n",
    "    eos.pressure = eos.pressure.subs(alpha, eos.alpha_function(reduced_temperature))\n",
    "\n",
    "    pressure, rho_V, rho_L = eos.calculate_pressure_and_densities(temperature)\n",
    "\n",
    "    rho_V_r = rho_V / rho_c\n",
    "    rho_L_r = rho_L / rho_c\n",
    "\n",
    "    gamma_VL_r = integral_VL = integrate_helholtz(eos, reduced_temperature, rho_V_r, rho_L_r, rho_inf=rho_V_r)\n",
    "\n",
    "    phi_values[reduced_temperature] = np.linspace(-phi_limits[reduced_temperature], phi_limits[reduced_temperature], num=25)\n",
    "\n",
    "    contact_angles[reduced_temperature] = []\n",
    "\n",
    "    for phi_r in phi_values[reduced_temperature]:\n",
    "        contact_angle = calculate_contact_angle(\n",
    "            phi_r, \n",
    "            gamma_VL_r=gamma_VL_r, \n",
    "            eos=eos, \n",
    "            reduced_temperature=reduced_temperature,\n",
    "            rho_V_r=rho_V_r,\n",
    "            rho_L_r=rho_L_r\n",
    "        )\n",
    "        contact_angles[reduced_temperature].append(contact_angle)\n",
    "\n",
    "    contact_angles[reduced_temperature] = np.array(contact_angles[reduced_temperature])\n",
    "\n",
    "    plt.plot(phi_values[reduced_temperature], \n",
    "             contact_angles[reduced_temperature], \n",
    "             marker='o', label='Reduced temperature = %.2f' %reduced_temperature)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.grid(True)\n",
    "\n",
    "ax.set_xlabel(xlabel='Specific Surface Energy [-]')\n",
    "ax.set_ylabel(ylabel='Contact Angle [deg]')\n",
    "ax.legend(loc='upper right')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Only Plot\n",
    "reduced_temperature_values = [0.8, 0.82, 0.87, 0.9]\n",
    "\n",
    "for reduced_temperature in reduced_temperature_values:\n",
    "    plt.plot(phi_values[reduced_temperature], \n",
    "             contact_angles[reduced_temperature], \n",
    "             marker='o', label='Reduced temperature = %.2f' %reduced_temperature)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.grid(True)\n",
    "\n",
    "ax.set_xlabel(xlabel='Specific Surface Energy [-]')\n",
    "ax.set_ylabel(ylabel='Contact Angle [deg]')\n",
    "ax.legend(loc='upper right')\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}