{
"cells": [
{
"cell_type": "markdown",
"id": "1afc4517",
"metadata": {},
"source": [
"# Thermodynamic Derivatives\n",
"\n",
"## Helmholtz energy derivatives\n",
"Thermodynamic derivatives are at the very heart of teqp. All models are defined in the form $\\alpha^r(T, \\rho, z)$, where $\\rho$ is the molar density, and z are mole fractions. There are exceptions for models for which the independent variables are in simulation units (Lennard-Jones and its ilk).\n",
"\n",
"Thereofore, to obtain the residual pressure, it is obtained as a derivative:\n",
"\n",
"$$ p^r = \\rho R T \\left( \\rho \\left(\\frac{\\partial \\alpha^r}{\\partial \\rho}\\right)_{T}\\right)$$\n",
"\n",
"and other residual thermodynamic properties are defined likewise.\n",
"\n",
"We can define the concise derivative\n",
"\n",
"$$ \\Lambda^r_{xy} = (1/T)^x(\\rho)^y\\left(\\frac{\\partial^{x+y}(\\alpha^r)}{\\partial (1/T)^x\\partial \\rho^y}\\right)$$\n",
"\n",
"so we can re-write the derivative above as \n",
"\n",
"$$ p^r = \\rho R T \\Lambda^r_{01}$$\n",
"\n",
"Similar definitions apply for all the other thermodynamic properties, with the tot superscript indicating it is the sum of the residual and ideal-gas (not included in teqp) contributions:\n",
"\n",
"$$\n",
"\\frac{p}{\\rho R T} = 1 + \\Lambda^r_{01}\n",
"$$\n",
"Internal energy ($u= a+Ts$):\n",
"$$\n",
"\t\\frac{u}{RT} = \\Lambda^{\\rm tot}_{10}\n",
"$$\n",
"Enthalpy ($h= u+p/\\rho$):\n",
"$$\n",
"\\frac{h}{RT} = 1+\\Lambda^r_{01} + \\Lambda^{\\rm tot}_{10}\n",
"$$\n",
"Entropy ($s\\equiv -(\\partial a/\\partial T)_v$):\n",
"$$\n",
"\\frac{s}{R} = \\Lambda^{\\rm tot}_{10}-\\Lambda^{\\rm tot}_{00}\n",
"$$\n",
"Gibbs energy ($g= h-Ts$):\n",
"$$\n",
"\t\\frac{g}{RT} = 1+\\Lambda^r_{01}+\\Lambda^{\\rm tot}_{00}\n",
"$$\n",
"Derivatives of pressure:\n",
"$$\n",
"\t\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_{T} = RT\\left(1+2\\Lambda^r_{01}+\\Lambda^r_{02}\\right)\n",
"$$\n",
"$$\n",
"\t\\left(\\frac{\\partial p}{\\partial T}\\right)_{\\rho} = R\\rho\\left(1+\\Lambda^r_{01}-\\Lambda^r_{11}\\right)\n",
"$$\n",
"Isochoric specific heat ($c_v\\equiv (\\partial u/\\partial T)_v$):\n",
"$$\n",
"\\frac{c_v}{R} = -\\Lambda^{\\rm tot}_{20}\n",
"$$\n",
"Isobaric specific heat ($c_p\\equiv (\\partial h/\\partial T)_p$; see Eq. 3.56 from Span for the derivation):\n",
"$$\n",
"\\frac{c_p}{R} = -\\Lambda^{\\rm tot}_{20}+\\frac{(1+\\Lambda^r_{01}-\\Lambda^r_{11})^2}{1+2\\Lambda^r_{01}+\\Lambda^r_{02}}\n",
"$$"
]
},
{
"cell_type": "raw",
"id": "4e621170",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
"In teqp, these derivatives are obtained from methods like \n",
"\n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_Arxy`\n",
"* :py:meth:`~teqp.teqp.AbstractModel.get_Ar06n`"
]
},
{
"cell_type": "markdown",
"id": "7622b9b7",
"metadata": {
"raw_mimetype": "text/restructuredtext"
},
"source": [
"where the ``A`` in this context indicates the variable $\\Lambda$ above. This naming is perhaps not ideal, since $A$ is sometimes the total Helmholtz energy, but it was a close visual mnemonic to the character $\\Lambda$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "2ced52fd",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:41:59.717164Z",
"iopub.status.busy": "2024-03-15T22:41:59.717003Z",
"iopub.status.idle": "2024-03-15T22:41:59.730446Z",
"shell.execute_reply": "2024-03-15T22:41:59.729997Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"'0.19.1'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import teqp\n",
"teqp.__version__"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ca35cc29",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:41:59.732347Z",
"iopub.status.busy": "2024-03-15T22:41:59.732058Z",
"iopub.status.idle": "2024-03-15T22:41:59.791168Z",
"shell.execute_reply": "2024-03-15T22:41:59.790693Z"
}
},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "8d254c40",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:41:59.794411Z",
"iopub.status.busy": "2024-03-15T22:41:59.793511Z",
"iopub.status.idle": "2024-03-15T22:41:59.797791Z",
"shell.execute_reply": "2024-03-15T22:41:59.797346Z"
}
},
"outputs": [],
"source": [
"Tc_K = [300]\n",
"pc_Pa = [4e6]\n",
"acentric = [0.01]\n",
"model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "49ee3453",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:41:59.800980Z",
"iopub.status.busy": "2024-03-15T22:41:59.800100Z",
"iopub.status.idle": "2024-03-15T22:41:59.805635Z",
"shell.execute_reply": "2024-03-15T22:41:59.805197Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"-0.06836660379313926"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"z = np.array([1.0])\n",
"model.get_Ar01(300,300,z)"
]
},
{
"cell_type": "markdown",
"id": "498a2350",
"metadata": {},
"source": [
"And there are additional methods to obtain all the derivatives up to a given order:\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "3ae755e1",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:41:59.808836Z",
"iopub.status.busy": "2024-03-15T22:41:59.807972Z",
"iopub.status.idle": "2024-03-15T22:41:59.813613Z",
"shell.execute_reply": "2024-03-15T22:41:59.813177Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([-6.96613834e-02, -6.83666038e-02, 2.53578225e-03, -1.57011622e-04,\n",
" 1.68186288e-05, -2.23059409e-06, 3.82592585e-07])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.get_Ar06n(300,300,z) # derivatives 00, 01, 02, ... 06"
]
},
{
"cell_type": "markdown",
"id": "cd644a3e",
"metadata": {},
"source": [
"But more derivatives are slower than fewer:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "4e0609ae",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:41:59.816848Z",
"iopub.status.busy": "2024-03-15T22:41:59.815971Z",
"iopub.status.idle": "2024-03-15T22:42:13.686778Z",
"shell.execute_reply": "2024-03-15T22:42:13.686315Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"601 ns ± 2.9 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.11 µs ± 2.43 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)\n"
]
}
],
"source": [
"%timeit model.get_Ar01(300,300,z)\n",
"%timeit model.get_Ar04n(300,300,z)"
]
},
{
"cell_type": "markdown",
"id": "d05e9070",
"metadata": {},
"source": [
"Note: calling overhead is usually on the order of 1 microsecond"
]
},
{
"cell_type": "markdown",
"id": "502e9830",
"metadata": {},
"source": [
"## Virial coefficients\n",
"\n",
"Virial coefficients represent the thermodynamics of the interaction of two-, three-, ... bodies interacting with each other. They can be obtained rigorously if the potential energy surface of interaction is fully known. In general, such a surface can only be constructed for small rigid molecules. Many simple thermodynamic models do a poor job of predicting the thermodynamics captured by the virial coefficients. \n",
"\n",
"The i-th virial coefficient is defined by \n",
"$$\n",
"\tB_i = \\frac{(\\alpha^r)^{(i-1)}}{(i-2)!}\n",
"$$\n",
"with the concise derivative term\n",
"$$\n",
"\t(\\alpha^r)^{(i)} = \\lim_{\\rho \\to 0} \\left(\\frac{\\partial^{i}\\alpha^r}{\\partial \\rho^{i}}\\right)_{T,\\vec{x}}\n",
"$$\n",
"\n",
"teqp supports the virial coefficient directly, there is the ``get_B2vir`` method for the second virial coefficient:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "720e120c",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:42:13.688912Z",
"iopub.status.busy": "2024-03-15T22:42:13.688583Z",
"iopub.status.idle": "2024-03-15T22:42:13.691881Z",
"shell.execute_reply": "2024-03-15T22:42:13.691436Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"-0.00023661263734465424"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.get_B2vir(300, z)"
]
},
{
"cell_type": "markdown",
"id": "2fba4369",
"metadata": {},
"source": [
"And the ``get_Bnvir`` method that allows for the calculation of higher virial coefficients:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "92d01992",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:42:13.693749Z",
"iopub.status.busy": "2024-03-15T22:42:13.693452Z",
"iopub.status.idle": "2024-03-15T22:42:13.696995Z",
"shell.execute_reply": "2024-03-15T22:42:13.696486Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{2: -0.0002366126373446542,\n",
" 3: 3.001768410777936e-08,\n",
" 4: -3.2409760373816364e-12,\n",
" 5: 3.961781646633723e-16,\n",
" 6: -4.5529239838367004e-20,\n",
" 7: 5.375927851118494e-24}"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.get_Bnvir(7, 300, z)"
]
},
{
"cell_type": "markdown",
"id": "d50998b1",
"metadata": {},
"source": [
"The ``get_Bnvir`` method was implemented because when doing automatic differentiation, all the intermediate derivatives are also retained.\n",
"\n",
"There is also a method to calculate temperature derivatives of a given virial coefficient"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "fa4ff386",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:42:13.698902Z",
"iopub.status.busy": "2024-03-15T22:42:13.698593Z",
"iopub.status.idle": "2024-03-15T22:42:13.701988Z",
"shell.execute_reply": "2024-03-15T22:42:13.701576Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"1.0095625628421257e-10"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.get_dmBnvirdTm(2, 3, 300, z) # third temperature derivative of the second virial coefficient"
]
},
{
"cell_type": "markdown",
"id": "c9a9823b",
"metadata": {},
"source": [
"## Isochoric Thermodynamics Derivatives\n",
"\n",
"In the isochoric thermodynamics formalism, the EOS is expressed in the Helmholtz energy density $\\Psi$ as a function of temperature and molar densities $\\vec\\rho$. This formalism is handy because it allows for a concise mathematical structure, well suited to implementation in teqp. For instance the pressure is obtained from (see https://doi.org/10.1002/aic.16074):\n",
"$$\n",
"p=-\\Psi + \\sum_{i=1}^N\\rho_i\\mu_i\n",
"$$\n",
"with the chemical potential $\\mu_i$ obtained from\n",
"$$\n",
"\\mu_i = \\left(\\frac{\\partial \\Psi}{\\partial \\rho_i}\\right)_{T,\\rho_{j\\neq i}}\n",
"$$\n",
"The molar densities $\\rho_i$ are related to the total density and the mole fractions:\n",
"$$\n",
"\\rho_i = x_i\\rho\n",
"$$\n",
"\n",
"In teqp, the isochoric derivative functions like ``get_fugacity_coefficients``, ``get_partial_molar_volumes`` take as arguments the temperature ``T`` and the vector of molar concentrations ``rhovec``=$\\vec\\rho$, which are obtained by multiplying the mole fractions by the total density.\n",
"\n",
"Example:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "d56640d3",
"metadata": {
"execution": {
"iopub.execute_input": "2024-03-15T22:42:13.703860Z",
"iopub.status.busy": "2024-03-15T22:42:13.703561Z",
"iopub.status.idle": "2024-03-15T22:42:13.748698Z",
"shell.execute_reply": "2024-03-15T22:42:13.748174Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([0.97884567, 0.99866748])"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"array([0.00470644, 0.00480351])"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"model = teqp.build_multifluid_model([\"CO2\",\"Argon\"], teqp.get_datapath())\n",
"T, rhovec = 300, np.array([0.3,0.4])*300 # K, mol/m^3\n",
"display(model.get_fugacity_coefficients(T, rhovec))\n",
"display(model.get_partial_molar_volumes(T, rhovec))"
]
}
],
"metadata": {
"celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}