Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add combination of multifluid and AC models (Wilson, COSMO-SAC) #155

Merged
merged 10 commits into from
Aug 29, 2024
254 changes: 254 additions & 0 deletions doc/source/models/MultifluidPlusAC.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,254 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "ef2685a7",
"metadata": {},
"source": [
"# Multifluid+Activity\n",
"\n",
"Following the derivations in [Jaeger et al.](https://doi.org/10.1016/j.fluid.2018.04.015), the multifluid plus activity coefficient model can be defined as a the conventional corresponding states portion\n",
"$$\n",
"\\alpha^{\\rm r}_{\\rm CS} = \\sum_i \\alpha^{\\rm r}_{oi}(\\tau, \\delta)\n",
"$$\n",
"plus the departure term coming from the activity coefficient model\n",
"$$\n",
"\\alpha^{\\rm Dep}(T,\\rho,\\bar x) = \\frac{\\ln(1+b_{\\rm mix}\\rho)}{\\ln(1+\\frac{1}{u})}\\left[\\frac{g^{\\rm E,R}_{\\rm GE}(T,\n",
"\\bar x)}{RT}-\\sum_{i=1}^Nx_i\\left(\\alpha^{\\rm r}_{oi}(\\delta_{\\rm ref},\\tau)-\\alpha^{\\rm r}_{oi}(\\delta_{i,ref}, \\tau_i)\\right)\\right]\n",
"$$\n",
"with \n",
"$$\n",
"b_{\\rm mix} = \\sum_i b_ix_i\n",
"$$\n",
"$$\n",
"\\frac{g^{\\rm E,R}_{\\rm GE}(T,\n",
"\\bar x)}{RT} = \\sum_ix_i\\ln\\left(\\gamma_i^{\\rm R}\\right)\n",
"$$\n",
"\n",
"$$\n",
"\\delta_{\\rm ref} = \\frac{1}{ub_{\\rm mix}\\rho_r(\\bar x)}\n",
"$$\n",
"$$\n",
"\\delta_{i,{\\rm ref}} = \\frac{1}{ub_i\\rho_{c,i}}\n",
"$$\n",
"and the conventional reducing function is used to define $T_{\\rm r}(\\bar x)$ and $\\rho_{\\rm r}(\\bar x)$.\n",
"\n",
"Any activity coefficient model can be used for $\\gamma_i^{\\rm R}$, but so far (as of version 0.22) only the Wilson model is available due to how well it works with the advanced cubic mixing rules."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "937c2216",
"metadata": {},
"outputs": [],
"source": [
"import io, teqp, numpy as np, CoolProp.CoolProp as CP\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e03c647",
"metadata": {},
"outputs": [],
"source": [
"# Four isotherms of experimental data for CO2 + Nitrogen from doi: 10.1016/j.fluid.2016.05.015\n",
"import pandas\n",
"dat = pandas.read_csv(io.StringIO(\"\"\"PointID y1 uy1 x1 ux1 p/bar up T/K\n",
"1 0.0274 0.0007 0.0068 0.0002 59.830 0.053 293.10 \n",
"2 0.0664 0.0014 0.0183 0.0004 64.864 0.080 293.10 \n",
"3 0.0978 0.0020 0.0298 0.0007 69.772 0.080 293.10 \n",
"4 0.1199 0.0024 0.0424 0.0009 74.737 0.080 293.10 \n",
"5 0.1219 0.0028 0.1132 0.0023 89.869 0.080 293.10 \n",
"6 0.1339 0.0024 0.0995 0.0022 89.198 0.080 293.10 \n",
"7 0.1399 0.0026 0.0943 0.0020 88.853 0.080 293.10 \n",
"8 0.1461 0.0027 0.0823 0.0019 86.962 0.080 293.10 \n",
"9 0.1466 0.0028 0.0778 0.0017 85.942 0.080 293.10 \n",
"10 0.1466 0.0028 0.0772 0.0016 85.868 0.080 293.10 \n",
"1 0.1378 0.0027 0.0159 0.0004 42.667 0.051 273.08 \n",
"2 0.2143 0.0038 0.0297 0.0007 49.547 0.051 273.08 \n",
"3 0.2612 0.0043 0.0411 0.0009 55.238 0.051 273.08 \n",
"4 0.3209 0.0049 0.0609 0.0013 65.069 0.088 273.08 \n",
"5 0.3554 0.0051 0.0786 0.0016 73.395 0.088 273.08 \n",
"6 0.3758 0.0052 0.0978 0.0019 81.061 0.088 273.08 \n",
"7 0.3903 0.0053 0.1190 0.0023 90.706 0.088 273.08 \n",
"8 0.3914 0.0053 0.1477 0.0028 100.966 0.088 273.08 \n",
"9 0.3879 0.0053 0.1614 0.0030 104.806 0.088 273.08 \n",
"10 0.3724 0.0052 0.1875 0.0033 110.846 0.088 273.08\n",
"11 0.3550 0.0051 0.2068 0.0036 114.105 0.088 273.08\n",
"12 0.2727 0.0044 0.2531 0.0041 118.020 0.088 273.08\n",
"13 0.3343 0.0049 0.2268 0.0038 116.295 0.088 273.08\n",
"1 0.2048 0.0038 0.0106 0.0003 25.754 0.050 253.05 \n",
"2 0.3019 0.0049 0.0217 0.0005 30.479 0.050 253.05 \n",
"3 0.4638 0.0056 0.0436 0.0010 45.352 0.050 253.05 \n",
"4 0.5319 0.0056 0.0647 0.0014 58.188 0.050 253.05 \n",
"5 0.5854 0.0054 0.1077 0.0021 78.315 0.084 253.05 \n",
"6 0.5979 0.0054 0.1497 0.0028 98.276 0.084 253.05\n",
"7 0.5898 0.0054 0.1801 0.0032 109.241 0.084 253.05\n",
"8 0.5042 0.0057 0.0570 0.0012 51.343 0.084 253.05\n",
"9 0.5644 0.0055 0.0861 0.0017 67.594 0.084 253.05\n",
"10 0.5949 0.0054 0.1267 0.0024 86.883 0.084 253.05\n",
"11 0.5826 0.0054 0.2015 0.0035 116.614 0.084 253.05\n",
"12 0.5537 0.0055 0.2431 0.0040 129.873 0.084 253.05\n",
"13 0.4973 0.0055 0.2971 0.0046 139.161 0.084 253.05\n",
"14 0.4971 0.0055 0.2972 0.0046 139.261 0.084 253.05\n",
"1 0.7076 0.0050 0.0257 0.0006 27.983 0.056 223.10\n",
"2 0.7774 0.0041 0.0522 0.0011 44.918 0.056 223.10\n",
"3 0.8077 0.0036 0.0930 0.0019 64.906 0.081 223.10\n",
"4 0.8131 0.0035 0.1261 0.0024 84.799 0.081 223.10\n",
"5 0.8057 0.0035 0.1584 0.0029 104.410 0.081 223.10\n",
"6 0.7843 0.0038 0.1982 0.0035 125.782 0.081 223.10\n",
"7 0.7533 0.0041 0.2380 0.0040 144.287 0.081 223.10\n",
"8 0.7150 0.0045 0.2813 0.0044 159.015 0.081 223.10\n",
"9 0.6942 0.0047 0.3064 0.0047 165.347 0.081 223.10\n",
"\"\"\"), sep='\\s+', engine='python')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "74404b8c",
"metadata": {},
"outputs": [],
"source": [
"Tc_K = np.array([304.21, 126.19])\n",
"pc_Pa = np.array([7.383e6, 3395800.0])\n",
"\n",
"def get_bi():\n",
" OmegaB = 0.08664 \n",
" R = 8.31446261815324\n",
" return (OmegaB*R*Tc_K/pc_Pa).tolist()\n",
"\n",
"donor = teqp.build_multifluid_model([\"CO2\", \"Nitrogen\"], teqp.get_datapath())\n",
"vc_m3mol = donor.get_vcvec()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ab83b1d7",
"metadata": {},
"outputs": [],
"source": [
"# Parameter set with linear mixing in the reducing function\n",
"BIP_linear = [{\n",
" \"Name1\": \"CarbonDioxide\",\n",
" \"Name2\": \"Nitrogen\",\n",
" \"betaT\": 1.0,\n",
" \"gammaT\": 0.5*(Tc_K[0]+Tc_K[1])/(Tc_K[0]*Tc_K[1])**0.5,\n",
" \"betaV\": 1.0,\n",
" \"gammaV\": 4*(vc_m3mol[0]+vc_m3mol[1])/(vc_m3mol[0]**(1/3) + vc_m3mol[1]**(1/3))**3,\n",
" \"F\": 0.0\n",
"}]\n",
"\n",
"# Parameter set with Lorentz-Berthelot\n",
"BIP_LorentzBerthelot = [{\n",
" \"Name1\": \"CarbonDioxide\",\n",
" \"Name2\": \"Nitrogen\",\n",
" \"betaT\": 1.0,\n",
" \"gammaT\": 1.0,\n",
" \"betaV\": 1.0,\n",
" \"gammaV\": 1.0,\n",
" \"F\": 0.0\n",
"}]\n",
"\n",
"def get_model(BIPset=None):\n",
" if BIPset is not None and BIPset == 'linear':\n",
" return teqp.make_model({\n",
" 'kind': 'multifluid',\n",
" \"model\":{\n",
" \"components\": [\"CO2\", \"Nitrogen\"],\n",
" \"root\": teqp.get_datapath(),\n",
" \"BIP\": BIP_linear\n",
" }})\n",
"\n",
" # The Wilson activity coefficient model is from Lasala et al.: https://doi.org/10.1016/j.fluid.2016.05.015\n",
" # and then combined with the \n",
" j = {\n",
" \"kind\": \"multifluid-activity\",\n",
" \"model\": {\n",
" \"multifluid\":{\n",
" \"components\": [\"CO2\", \"Nitrogen\"],\n",
" \"root\": teqp.get_datapath(),\n",
" \"BIP\": BIPset\n",
" },\n",
" \"activity\": {\n",
" \"aresmodel\": {\n",
" \"type\": \"Wilson\", \n",
" \"m\": [[0.0, -3.4768], [3.5332, 0.0]], \n",
" \"n\": [[0.0, 825], [-585, 0.0]], \n",
" \"b\": get_bi()\n",
" },\n",
" \"options\": {\"u\": 1.17, \"b\": get_bi()}\n",
" }\n",
" }\n",
" }\n",
"\n",
" return teqp.make_model(j) if BIPset is not None else donor"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "500aa1ae",
"metadata": {},
"outputs": [],
"source": [
"labels = ['linear+Wilson','Lorentz-Berthelot+Wilson','GERG-2008','linear']\n",
"for imod, model in enumerate([get_model(BIP_linear), get_model(BIP_LorentzBerthelot), get_model(), get_model('linear')]):\n",
" lw = 0.5+imod\n",
" for T in [223.15, 253.05, 273.08, 293.1]:\n",
" ipure = 0\n",
"\n",
" [rhoL0, rhoV0] = [CP.PropsSI('Dmolar','T',T,'Q',Q,'CO2') for Q in [0,1]]\n",
"\n",
" rhovecL0 = np.array([0.0, 0.0]); rhovecL0[ipure] = rhoL0\n",
" rhovecV0 = np.array([0.0, 0.0]); rhovecV0[ipure] = rhoV0\n",
"\n",
" J = model.trace_VLE_isotherm_binary(T, rhovecL0, rhovecV0)\n",
" df = pandas.DataFrame(J)\n",
" plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6,'k', lw=lw, label=labels[imod] if T == 223.15 else '')\n",
" plt.plot(df['xV_0 / mole frac.'], df['pV / Pa']/1e6,'k', lw=lw)\n",
" \n",
"plt.title('CO$_2$ (1) + N$_2$ (2)')\n",
"plt.plot(1-dat['x1'], dat['p/bar']/10, 'o')\n",
"plt.plot(1-dat['y1'], dat['p/bar']/10, '^')\n",
"plt.legend(loc='best')\n",
"plt.gca().set(xlabel='$x_1$ / mole frac.', ylabel='$p$ / MPa', ylim=(0, 25))\n",
"plt.savefig(\"multifluid_Wilson_CO2_N2.pdf\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e1b7813a",
"metadata": {},
"source": [
"so we can see that the activity coefficient models are not too bad, clearly an improvement over pure linear mixing, but the model from GERG-2008 is unsurprisingly the most accurate"
]
}
],
"metadata": {
"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.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
1 change: 1 addition & 0 deletions doc/source/models/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -28,5 +28,6 @@ Models
GERG
ECS
MultifluidPlusAssociation
MultifluidPlusAC

IdealGas
Loading
Loading