Newer
Older
"outputs": [
{
"data": {
"text/plain": [
"'0.9.4.dev0'"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import teqp\n",
"teqp.__version__"
]
},
{
"cell_type": "markdown",
"id": "8218498b",
"metadata": {},
"source": [
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
"# Critical curves & points\n",
"\n",
"\n",
"## Pure Fluids\n",
"\n",
"Solving for the critical point involves finding the temperature and density that make\n",
"$$\n",
"\\left(\\frac{\\partial p}{\\partial \\rho}\\right)_T = \\left(\\frac{\\partial^2 p}{\\partial \\rho^2}\\right)_T = 0\n",
"$$\n",
"by 2D non-linear rootfinding. Newton steps are taken, and the analytic Jacobian is used (thanks to the ability to do derivatives with automatic differentiation). This is all handily wrapped up in the\n",
"``solve_pure_critical`` method which requires the user to provide guess values for temperature and density"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "46657a96",
"metadata": {},
"outputs": [],
"source": [
"# Values taken from http://dx.doi.org/10.6028/jres.121.011\n",
"modelPR = teqp.canonical_PR([190.564], [4599200], [0.011])\n",
"\n",
"# Solve for the critical point from a point close to the critical point\n",
"T0 = 192.0\n",
"# Critical compressibility factor of P-R is 0.307401308698.. (see https://doi.org/10.1021/acs.iecr.1c00847)\n",
"rhoc = (4599200/(8.31446261815324*190.564))/0.3074\n",
"rho0 = rhoc*1.2345 # Perturb to make sure we are doing something in the solver\n",
"teqp.solve_pure_critical(modelPR, T0, rho0)"
]
},
{
"cell_type": "markdown",
"id": "e15eeded",
"metadata": {},
"source": [
"\n",
"A pure fluid has a single vapor-liquid critical point, but mixtures are different:\n",
"\n",
"* They may have multiple (or zero!) critical points for a given mixture composition\n",
"* The critical curves may not emanate from the pure fluid endpoints\n",
"\n",
"When it comes to critical points, intuition from pure fluids is not helpful, or sometimes even counter-productive. \n",
"\n",
"teqp has methods for working with the critical loci of binary mixtures (only binary mixtures, for now) and especially, methods for tracing the critical curves emanating from the pure fluid endpoints.\n",
"\n",
"The tracing method in teqp is based explicitly on the isochoric thermodynamics formalism introduced by Ulrich Deiters and Sergio Quinones-Cisneros. It uses the Helmholtz energy density as the fundamental potential and all other properties are derived from it. For critical curves it is based upon the integration of sets of ordinary differential equations; the differential equations are in the form of derivatives of the molar concentrations of each component in the mixture with respect to an integration variable. The set of ODE is then integrated.\n",
"\n",
"Here is an example of the construction of the critical curves emanating from the pure fluid endpoints for the mixture nitrogen + ethane."
]
},
{
"cell_type": "code",
"id": "81619ae2",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgYAAAG8CAYAAAClhm0uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4zElEQVR4nO3dd3gc5b328e9PvVqyJBfZlhtyxTa2kU2vIWA6hBKMIbRASELqCYG095DknJxDck5OQggBBwi9hYRiQugYMMW9927LTZZsSVaXdp/3j5U3wkjCZaXZnb0/16VrV7PS+vZ4rLn1zMwz5pxDREREBCDB6wAiIiISPVQMREREJEzFQERERMJUDERERCRMxUBERETCVAxEREQkTMVAREREwlQMREREJCzqi4GZjTKz+83seTP7utd5RERE/MyTYmBmD5tZmZktO2D5FDNbbWbrzOxOAOfcSufcrcCVQIkXeUVEROKFVyMGjwBT2i4ws0Tgj8C5wGhgqpmNbn3tImAW8Hb3xhQREYkvnhQD59z7wJ4DFk8G1jnnNjjnmoBngItbv/5l59yJwLTuTSoiIhJfkrwO0EZ/YGubz0uB48zsdOBLQCrwakffbGa3ALcAZGZmHjty5MguCyoi3a++KcC+xhZqWz8cYEBGShJZaUlkpSaRnpKIeR1UxCPz588vd871OtL3iaZi0N7/Z+ecmwnM/Lxvds5NB6YDlJSUuHnz5kU0nIh0H+ccmyvq+GBdOR+uLeej9eXUNrSQAEwu7MHJxfmcVFzA5CF5ZKRE048xEe+Y2eZIvE80/Y8qBYrafD4A2O5RFhHpZuU1jXy0voIP15Yza1052yrrAeiXk8aUMX05qbiAE48qoFd2qsdJRfwtmorBXGCYmQ0BtgFXAVd7G0lEukpdUwtzNu7hw3XlzFpXwcod1QD0SEvixKMKuPX0ozi5uIDB+RmY6QCBSHfxpBiY2dPA6UCBmZUC/+6ce8jMbgNeBxKBh51zy73IJyJdY2N5LW+v3MXbK8uYv3kvTYEgKYkJHDuoJ7efM4KTigsY2z+HxAQVARGveFIMnHNTO1j+Kp2cYCgisaUlEGTe5r3hMrChvBaAEX2yuf6kwZxcXMCkwXmkpyR6nFRE9oumQwlHzMwuBC4sLi72OopI3Kqqa2bmmjLeWVXGzNW7qapvJjnROH5oPtedOJgzR/amKC/D65gi0gFfFQPn3AxgRklJyc1eZxGJJxt21/DOqjLeWrmLuZv2Egg68jNT+OLoPnxhZG9OGd6LrFRf/bgR8S39TxWRQ9bZIYKvnTqUL4zqw/iiXJ0rIBKDVAxE5KDsP0Tw9soyZq4uo7qhhZTEBI4bmqdDBCI+omIgIh3asLuGt1eGDhHM2/yvQwRnH92Xs0b15uRhOkQg4jf6Hy0iYS2BIHM37eWdVZ8+RDCybza3njaUM0fqEIGI3/mqGOiqBJFD19EhguOPyuf6kwZzxggdIhCJJ+ac8zpDxOleCSKd21vbxGvLdzJj8XZmb9xDIOgoyErhjBG9+YIOEYjEJDOb75wrOdL30f98kThR09jCWyt28fLi7by/ZjctQcfQgkxuPa31KoIBuSToEIFI3FMxEPGxhuYAM1fvZsbi7by9ahcNzUH65aRx08lDuPCYfhzdr4fuQyAin6JiIOIzLYEgH66v4OVF23lj+U72NbaQn5nClSVFXHRMPyYO7KmRARHpkIqBiA8Eg455m/fy8uJtvLp0J3tqm8hOS2LKmL5cNL4fJwzNJykxweuYIhIDVAxEYpRzjmXbqnl58TZeWbKDHVUNpCUncNaoPlx4TD9OH9GL1CTdnEhEDo2KgUiMWVe2j5cXbWfGkh1sLK8lOdE4bXgv7jx3JGeN6kOmriYQkSPgq58gmsdA/GrrnjpmLNnOjMU7WLmjmgSDE47K52unDmXKmL7kZqR4HVFEfELzGIhEqbJ9DfxjyQ5mLN7Ogi2VAEwcmMuFx/Tj/HGF9M5O8zagiEQVzWMg4kOBoOPdVWU8OXsz763ZTdDBqMIe/HDKCC4c108zEIpIl1MxEIkCZfsaeG7uVp6es5VtlfX0zk7lG6cXc/H4fgzrk+11PBGJIyoGIh5xzvHJhj08MXszry/bSUvQcVJxPj89fxRnje5Dsi4vFBEPqBiIdLOq+mb+Nr+UJ2dvZv3uWnLSk7nuxMFMO24gQ3tleR1PROKcioFIN1lSWskTn2zm5cXbaWgOMr4ol/+54hguGFdIWrLmGxCR6KBiINKF6psCzFi8nSdmb2ZJaRXpyYlcOqE/044bxJj+OV7HExH5DF8VA81jINFiXdk+nvhkC39bUMq+hhaG9c7i5xcdzaUT+9MjLdnreCIiHfJVMXDOzQBmlJSU3Ox1Fok/TS1B3lixkyc+2cwnG/aQnGicO6aQa44fxKTBPXUXQxGJCb4qBiJeKN1bxzNztvLM3K2U1zQyoGc6P5wygitLiijISvU6nojIIVExEDkMgaDj/TW7eeKTzby7ugwHnDmiN9ccP4hTh/ciUbc1FpEYpWIgcgjKaxp5du5Wnp6zhdK99RRkhSYiumpyEQN6alZCEYl9KgYin8M5x9xNe3n8k828tmwHzQHH8UPzuPPckZw9ui8pSZqISET8Q8VApAPOOd5bs5vfv72WhVsqyU5LYtpxg7jm+IEU99Y0xSLiTyoGIgdwzjFzdagQLNpaSf/cdH55yRgunziA9BRNRCQi/qZiINLKOcc7q8r4/dtrWVJaRf/cdH516VguP3aADheISNxQMZC455zjrZVl3PP2WpZuq6IoL527LxvLpRNUCEQk/qgYSNwKBh1vrNjFPW+vZcWOagblZ/Dry8dx6YT+urOhiMQtXxUDTYksByMYdLy+fCe/f3stq3buY3B+Bv9zxTFcMr4fSSoEIhLnzDnndYaIKykpcfPmzfM6hkSZYNDxz2U7+cM7oUIwtCCT284s5qJjVAhEJPaZ2XznXMmRvo+vRgxE2hMIOl5duoM/vLOWNbtqOKpXJr+/ajwXjOunGQpFRA6gYiC+FQg6XlmynT+8s451ZTUM653FPVMncP7YQhUCEZEOqBiI77QEgryyJDRCsH53LcP7ZHHv1RM4b0whCSoEIiKdUjEQ32gJBHlp0XbufXcdG8trGdk3m/umTWTK0X1VCEREDpKKgfjCR+vL+dmLy1i/u5ZRhT24/5qJnD1ahUBE5FCpGEhM272vkf/8xwpeXLSdgXkZ3H/NsZxzdB/MVAhERA6HioHEpEDQ8dScLfz6tVU0Ngf59pnFfOOMYtKSdS8DEZEjoWIgMWfZtip+8sJSFpdWcVJxPr+4eAxH9cryOpaIiC+oGEjMqG5o5rdvrOGxjzeRl5nK768az0XH9NNhAxGRCFIxkKjnnGPGkh388pUVlNc08pXjB/H9s0eQk57sdTQREd/xVTHQvRL8Z2N5LT97cRmz1pUztn8OD11XwrgBuV7HEhHxLd0rQaJSQ3OA+2au5/6Z60lNSuD2KSOYdtwgzVgoItIB3StBfOu9Nbv5fy8tY3NFHReP78dPzh9F7+w0r2OJiMQFFQOJGjurGvjlKyv4x9IdDC3I5MmvHsdJxQVexxIRiSsqBuK5lkCQxz7ezG/fXENzIMi/fXE4t5w2lNQkzUkgItLdVAzEUwu37OUnLyxjxY5qTh/Ri19cNIaB+RlexxIRiVsqBuKJqrpm7n59FU/P2UKf7DT+NG0iU8b01ZwEIiIeUzGQbuWc44WF2/jPf6yksr6Zm04awne/OJysVG2KIiLRQD+NpdvUNbXw478v5cVF25k4MJcnLh3LqMIeXscSEZE2VAykW6zfXcPXn5jP2rIafnD2cL5xerFuiSwiEoVUDKTLvbp0Bz98fgkpSQk8fuNxnDxMlyCKiEQrFQPpMs2BIHf/cxUPztrIhIG53DdtIoU56V7HEhGRTqgYSJfYVd3AbU8tYO6mvVx/4mB+fN4oUpISvI4lIiKfQ8VAIu6TDRXc9tRC6ppauGfqBC46pp/XkURE5CCpGEjEOOd44P0N/Ob11QzOz+Dpm49jWJ9sr2OJiMghUDGQiKhuaOYHzy3mjRW7OH9cIXdfNk5zE4iIxCBf/eQ2swuBC4uLi72OEldW7qjm60/Mp3RvPf/vgtHccNJgzWAoIhKjfHU2mHNuhnPulpycHK+jxI2/zS/l0vs+pL45wDO3HM+NJw9RKRARiWG+GjGQ7tPQHODnM1bw9JwtnDA0nz9cPYGCrFSvY4mIyBFSMZBDVravga8+Oo8lpVV84/Sj+P4Xh5OU6KvBJxGRuKViIIdke2U90x6cza7qBv78lRK+OLqP15FERCSCVAzkoG2uqOXqP8+muqGZx2+azLGD8ryOJCIiEaZiIAdl7a59THtwNs2BIE/ffDxj+usETxERP1IxkM+1fHsV1z40h8QE49mvncBwTVokIuJbOmNMOrVgy16mTv+E9ORE/qpSICLiexoxkA59vL6Crz46l17ZqTzx1eMY0DPD60giItLFVAykXe+uLuPWx+czMC+DJ796HL17pHkdSUREuoGKgXzGa8t28K2nFzK8TzaP33QceZkpXkcSEZFuomIgn/LCwlJ+8NclHDMgh7/cMJmc9GSvI4mISDdSMZCwFxaW8v3nFnP8kHwevK6ETN0dUUQk7ugnvwAwd9Mefvj8Ek4Yms/D108iLTnR60giIuIBXa4obN1Tx9cen09Rzwz+NO1YlQIRkTimYhDn9jU0c9OjcwkEHQ9dP4mcDJ1TICISz1QM4lgg6Pj20wvZsLuWP02byJCCTK8jiYiIx3SOQRz7r1dX8u7q3fznpWM4sbjA6zgiIhIFNGIQp56Zs4UHZ23k+hMHM+24QV7HERGRKKFiEIc+Xl/BT19cxqnDe/HT80d5HUdERKKIikGc2VxRy9efnM/ggkzuvXoCSYnaBERE5F+0V4gj1Q3N3PjIXAAeuq6EHmm6AkFERD7NV8XAzC40s+lVVVVeR4k6LYEgtz21kM0Vddx/zbEMytcVCCIi8lm+KgbOuRnOuVtycnK8jhJ1fvP6at5fE7oC4fih+V7HERGRKOWrYiDtW7BlL9M/2MDUyQP58qSBXscREZEopmLgc40tAe54fgmFPdL48XkjvY4jIiJRThMc+dwf313P2rIa/nLDJLJ1sqGIiHwOjRj42Mod1dz37jq+NKE/Z4zo7XUcERGJASoGPtUSCHLH35aQk57Mzy4Y7XUcERGJETqU4FMPzdrIktIq7r16Aj0zU7yOIyIiMUIjBj60sbyW3765hrNH9+H8sYVexxERkRiiYuAzwaDjjr8tISUpgV9eMgYz8zqSiIjEEBUDn3lqzhbmbNzDT88fRZ8eaV7HERGRGKNi4COVdU3c/c9VnFScz5UlRV7HERGRGKRi4CMPf7iJfY0t/OyC0TqEICIih0XFwCeqG5p55MONnHN0H0b27eF1HBERiVEqBj7x+MebqW5o4VtnDvM6ioiIxDAVAx+oa2rhwQ82cMaIXozprztLiojI4VMx8IEnP9nC3rpmbtNogYiIHCEVgxjX0Bxg+gcbOKk4n2MH9fQ6joiIxDgVgxj37Nyt7N7XyG1naLRARESOnIpBDGtqCXL/e+uZNLgnxw/N8zqOiIj4gIpBDPv7glJ2VDVw25nDNG+BiIhEhIpBjGoJBLlv5nqOGZDDqcMKvI4jIiI+oWIQo95csYste+r45hnFGi0QEZGIUTGIUTOWbKcgK5UvjOrjdRQREfERFYMYVNvYwjuryjhvbF8SEzRaICIikaNiEIPeXlVGQ3OQC8b18zqKiIj4jIpBDHpl8Xb69EilRBMaiYhIhKkYxJh9Dc3MXLOb88YWkqDDCCIiEmEqBjHmzRW7aGrRYQQREekaKgYx5pUlO+ifm87EgbleRxERER9SMYghVXXNfLB2N+ePK9TcBSIi0iVUDGLI68t30hxwXDCu0OsoIiLiUyoGMeT15TspyktnbP8cr6OIiIhPqRjECOcc87fs5cShBTqMICIiXUbFIEZsqqijsq6ZCTrpUEREulBMFAMzu8TM/mxmL5nZ2V7n8cLCLXsBGK9iICIiXcizYmBmD5tZmZktO2D5FDNbbWbrzOxOAOfci865m4HrgS97ENdzi7ZWkpmSyLDe2V5HERERH/NyxOARYErbBWaWCPwROBcYDUw1s9FtvuSnra/HnYVbKhk3IFc3TRIRkS7lWTFwzr0P7Dlg8WRgnXNug3OuCXgGuNhC7gb+6Zxb0N1ZvdbQHGDljmqdXyAiIl0u2s4x6A9sbfN5aeuybwFnAZeb2a3tfaOZ3WJm88xs3u7du7s+aTdatq2KlqBjfFGu11FERMTnkrwOcID2xsmdc+4e4J7OvtE5Nx2YDlBSUuK6IJtnFm2tBHTioYiIdL1oGzEoBYrafD4A2O5RlqixcEsl/XPT6Z2d5nUUERHxuWgrBnOBYWY2xMxSgKuAlz3O5LkVO6oZN0CzHYqISNfz8nLFp4GPgRFmVmpmNznnWoDbgNeBlcBzzrnlXmWMBs45tlfWM6BnutdRREQkDnh2joFzbmoHy18FXu3mOFGrqr6ZxpYgfXroMIKIiHS9aDuUcETM7EIzm15VVeV1lIjZWd0AoGIgIiLdwlfFwDk3wzl3S06Of47H76wKFYO+OSoGIiLS9XxVDPyorLoRgL4aMRARkW6gYhDl9h9K6N0j1eMkIiISD1QMotzO6gbyMlNITUr0OoqIiMQBFYMot6uqQSceiohIt1ExiHI7qxvoq8MIIiLSTXxVDPx4ueKu6gZdkSAiIt3GV8XAb5crOucor2miIEsjBiIi0j18VQz8Jth6j8jkRP0ziYhI99AeJ4oFXagZJLR3M2oREZEuoGIQxfYXAzM1AxER6R4qBlGstReQoGIgIiLdRMUgiulQgoiIdDdfFQO/Xa4Y1IiBiIh0M18VA79drvivcww8DiIiInHDV8XAb1ww9KgRAxER6S4qBlHMoREDERHpXioGUSwlKfTP09gS9DiJiIjECxWDKJaenEhKUgJ765q8jiIiInFCxSCKmRm56clU1TV7HUVEROKEikGUy81IplLFQEREuomKQZTLTU/RoQQREek2vioGfpvgCEIjBlX1GjEQEZHu4ati4LcJjkCHEkREpHv5qhj4UW5GCpX1OpQgIiLdQ8UgyuVmJNPQHKShOeB1FBERiQMqBlEuLyMFgPKaRo+TiIhIPFAxiHJFeRkAbNlT53ESERGJByoGUW5QfqgYbK5QMRARka6nYhDlCnPSSUlMUDEQEZFuoWIQ5RITjAF56WyuqPU6ioiIxAEVgxgwOD+TTRoxEBGRbuCrYuDHmQ8hdJ7B5opanHNeRxEREZ/zVTHw48yHEBoxqGsKUF6jiY5ERKRr+aoY+NX+KxM26TwDERHpYioGMWB4n2wAVu6o9jiJiIj4nYpBDCjMSaMgK5XFW/117oSIiEQfFYMYYGaMG5DD0m2VXkcRERGfUzGIEWP757CurIbaxhavo4iIiI+pGMSIY4pyCDpYvl3nGYiISNdRMYgRY/vnArCktNLTHCIi4m9Jh/NNZtYTGAak7V/mnHs/UqHks3plp9IvJ40lpToBUUREus4hFwMz+yrwHWAAsAg4HvgYODOiyeQzxg7I0YiBiIh0qcM5lPAdYBKw2Tl3BjAB2B3RVNKucQNy2VRRR1Vds9dRRETEpw6nGDQ45xoAzCzVObcKGBHZWNKecQNCUz0v0WWLIiLSRQ6nGJSaWS7wIvCmmb0EbI9kqMPl15so7TdhYE+SEoyP1ld4HUVERHzqkIuBc+5S51ylc+4u4GfAQ8AlEc51WPx6E6X9slKTmDAwl1lry72OIiIiPnXQxcDM0szsu2Z2r5l9zcySnHPvOededs7ptn/d5OTiXizbXsXeWq1yERGJvEMZMXgUKAGWAucC/9sliaRTJw8rwDn4cL1GDUREJPIOpRiMds5d45x7ALgcOKWLMkknjhmQQ3Zakg4niIhIlziUYhC+Rs45pwn7PZKUmMAJQ/P5YG05zjmv44iIiM8cSjE4xsyqzWyfme0Dxu1/bmaawL8bnTKsgG2V9WyqqPM6ioiI+MxBz3zonEvsyiBy8E4e1guAWWt3M6Qg0+M0IiLiJwddDMzs5c5ed85ddORx5GAMzs+gf246768t59oTBnsdR0REfORQ7pVwAlAKPAXMBqxLEsnnMjO+MKo3z83bSl1TCxkph3UvLBERkc84lHMM+gI/AsYAvwe+CJS3zmXwXleEk46dO6aQhuYgM1frNhUiIhI5B10MnHMB59xrzrnrCN1RcR0w08y+1WXppEOTh+SRn5nCq0t3eB1FRER85JDGoM0sFTgfmAoMBu4B/h75WPJ5EhOMc8b05cWF22hoDpCWrHNDRUTkyB3KlMiPAh8BE4GfO+cmOed+6Zzb1mXppFPnjSmkrinAe2t0OEFERCLjUM4xuBYYDnwH+Kh1ToNqzWPgneOG5tEzI5l/6nCCiIhEyKHMY3A4t2iWLpScmMDZo/vyj6U7aGwJkJqkwwkiInJkfLWzN7MLzWx6VVWV11G6zXnjCqlpbOHtlWVeRxERER/wVTFwzs1wzt2Sk5PjdZRuc3JxAf1z03lq9havo4iIiA/4qhjEo8QE48uTipi1rpzNFbVexxERkRh3KFclnGBmmu0wCl1ZUkRigvH0nK1eRxERkRh3KCMG1wHzzewZM7vezPp2VSg5NH1z0jhzZG+en7+Vppag13FERCSGHcrMh7c65yYCdwE9gUfM7GMz+5WZnWpmOiXeQ1cfN5DymibeWLHT6ygiIhLDDvkcA+fcKufc/znnpgBnArOAKwjdWEk8cuqwXjoJUUREjtgRnXzonKt3zr3qnPuWc64kUqHk0CUmGFMnF/HR+go2luskRBEROTy6KsFHriwpIinBeGaORg1EROTwqBj4SO8eaZw1qg9/nV9KY0vA6zgiIhKDDrkYmNm5ZjbbzFab2XNmdkJXBJPDc/VxA9lT28Try3d5HUVERGLQ4YwY3Ad8HzgemA78xsymRjSVHLaTiwsYmJfBIx9uxDnndRwREYkxh1MMdjnnPnTO7XXOvQWcA/wkwrnkMCUkGDefOpQFWyr5YG2513FERCTGHE4x2GRm/2FmKa2fNwP7IphJjtCVJQPol5PG795ao1EDERE5JIdTDBzwJWCrmc0C1gEzzWxYRJPJYUtNSuSbZxazYEsl72vUQEREDsHhTHA01Tk3GhgEfBf4OZAJPGhmmqw/SlxxbBH9c9P5vzc1aiAiIgfvsC9XdM41OOfmOececs592zl3mnOuKJLh5PClJCVw25nFLNpaycw1u72OIyIiMULzGPjY5ccOYEDPdH6nUQMRETlIKgY+lpyYwLfOLGZxaRXvri7zOo6IiMQAFQOf+9LEAQzMy+B3b63VqIGIiHwuFQOfS04MnWuwpLSKN1ZoNkQREemcikEc+NKE/gzrncUvZqygtrHF6zgiIhLFfFUMzOxCM5teVVXldZSokpSYwK++NJZtlfX835trvI4jIiJRzFfFwDk3wzl3S05OjtdRos6kwXlMnTyQhz/cyLJtKk4iItI+XxUD6dydU0aSl5nKj/6+lJZA0Os4IiIShVQM4khORjL/fuFolm6r4rGPN3sdR0REopCKQZy5YFwhp4/oxf++sZrtlfVexxERkSijYhBnzIxfXjyGgHP8v5eWa24DERH5FBWDOFSUl8H3zhrOWyt38fz8Uq/jiIhIFFExiFM3nTyEE4bm89MXl+kqBRERCVMxiFNJiQn84eoJ5GWmcOsT89lb2+R1JBERiQIqBnGsICuV+6ZNpKy6ke88u4hAUOcbiIjEOxWDODdhYE/uuuho3l+zW7MiioiIioHA1MlFXFkygHvfXcfry3d6HUdERDykYiCYGb+4eAzjBuRw21MLmLF4u9eRRETEIyoGAkBaciKP33gcE4p68q2nF/LwrI1eRxIREQ+oGEhYTkYyj900mXOO7sMvXlnBf/1zJUGdkCgiEldUDORT0pITuW/asVxz/EAeeG8D//bXxTS16IZLIiLxIsnrABJ9EhNC0yb37ZHG/7yxhvKaRv50zbFkpWpzERHxO/2kl3aZGbedOYxe2an8+IVlTJ3+CQ9fP4le2aleR4s7gaCjtqmFmoYWahpDH80tQQLOEQxCwDkCwSCBYOhrg87REnQEg45A0LW+HlqelGAkJiSQlGAkJVr485SkBNKTE0lL3v+YSHpK62NyIokJ5vVqEJFuomIgnfrypIEUZKXyzacWcPn9H/HYjZMZlJ/pdayYEwg69tQ2Ubavgd37Ginb10h5TSPV9S3UNv5rh1/T0BIuAfsaQ6/VNQW8jk9qUgLZacn0SE+iR1oyPdKT6ZGWRI/0ZLLTQsvyMlMoyEolPyuFgszQY6ZGmURijvnx7nolJSVu3rx5XsfwlQVb9nLTI3NJMOOOc0dy2cQB+i0SqG8KtO7o/7XDb+/zitqmdmeWTE40slKTyEpLIjMliey0JDJTk8hKbX2eEnotq3XZ/tdSkhJIMCMxwUhMgMSEBBLNSEiApIQEEhNo83roI8EsNIIQdDQHggSCoZGFQNDR2BKgoTlIfVOA+uYADa0foedBahtbqG5oobqhmer6ZqobWtjX0Ex1fQvV9c00Bdo/DyU9OZH8rBTys1LpnZ1Kv5w0CnPTKcxJo39uOoW56fTJTiUpUac7iRwpM5vvnCs54vdRMZCDtX53Dd97dhFLSqsY1juLH5wzgrNH98HMvwWhoTlA6d46tu6tp3RP6HHrnjpK99azdW8dlXXNn/meBAtNN927Ryq9slLpnZ1Gr+zQ572zU0PPs9MoyEolPSXRg79V5DU0B9hT20R5TSMVNaHH8pomKmpCpai8ppFd1Q3sqGxgX2PLp743waB3dhqFuaGyMDg/k8EFmQzOz2BwQSb5mSm+3sZEIkXFoBMqBl3HOcdry3bymzdWs2F3LeOLcrljykhOOCrf62iHraaxhU3ltWyuqGNTRS2bymtDjxV17N7X+KmvTUlKYEDPdAb0zKCoZzr9ctM/tbPvlZ1KXmaKRlM6Ud3QzI7KBrZX1bOjsoEdVfVsb33cVllP6d76T42uZKcmhYpCQSZDWstCce8shvXO9k2xEokEFYNOqBh0vZZAkOfnl/K7t9ays7qBU4f34vazRzB2QI7X0dq1f+e/qSJUADaW17K5opaN5XWU13x65987OzX8G+vAvIxQCchLp6hnBgVZqSRop9+lmgNBtu6pC/87baqoDT9u21vP/s5gBoPzMxnRJ5sRfbMZ2Tf0OCg/U8VM4pKKQSdUDLpPQ3OAxz7exB/fXU9VfTP9c9OZPCSPyUPymDQ4j6N6ZXbbMHDbnX/osS78eODOv0+PVAblZzIkP5NBBRmhx/xMBhdkkJGiE+aiVVNLkC176lhXto9VO/exemfocVNFLft/lKUlJzCsdzajC3swfmAu44tyGd4nW2VBfE/FoBMqBt2vqr6ZFxduY/bGCuZs3EN5TRMA+Zkp4ZIweUgeowp7HPYP6OZAkL11TZRVN7az868N/5n7td35tz1mPShfO3+/qW8KsLZNWVi9cx/LtleFzwHJTElk7IAcxhf1ZMLAXCYU5dK7R5rHqUUiS8WgEyoG3nLOsbG8ljkb94Q+Nu2hdG89EDpeXNwni7Sk0DXzqQc+JicSDDr21DWxt7aJvXXN7K1rYk9tE/saWj7zZ/XpkRo6WU07fzmAc45NFXUs2rqXRVsqWbi1khXbq2lpPRbRLyeNCQN7Mr4ol/EDcxnbP4e0ZJ2zILFLxaATKgbRZ3tlPXM37WH2xj1srqilsTlIY0swfJlcY0uAxpYgDc0BEszomZFCXmYKuRmh6+N7ZqS0LkumV3ZoJEA7fzlUDc0Blm+vYuGWShZtrWThlkq2VYZKa0pSAiWDenLysAJOLi7g6H45OvwgMUXFoBMqBiJysMr2NbBoSyWzN+7hw3XlrNq5D4DcjGROPCqfk4t7ccqwAoryMjxOKtK5SBUD/bolInGtd3YaZx/dl7OP7guEisJH6yqYta6cWWvLeXXpTgAG5mWERxNOPCqf3IwUL2OLdBmNGIiIdMA5x/rdtXy4rpwP1pbzyYYKahpbMIMJRbmcO6aQKWP6ajRBooIOJXRCxUBEukJLIMji0ko+WFvOWyt3sWxbNQBj+vfg3DGFnDumL0N7ZXmcUuKVikEnVAxEpDtsqajjteU7eHXpThZtrQRgRJ9szh3bl3PHFDK8T5amc5Zuo2LQCRUDEelu2yvreW3ZTl5btpO5m/fgHAwtyGTKmL6cN7aQo/v1UEmQLqVi0AkVAxHxUtm+Bl5fvovXlu3gkw17CAQdowp7MHVyERcf05+cjGSvI4oPqRh0QsVARKLFntom/rFkO8/M3cry7dWkJiVw3thCrppUxOQheRpFkIhRMeiEioGIRKNl26p4Zu4WXlq4nX2NLQwpyOTLk4q4bOIAemWneh1PYlzcFAMzGwr8BMhxzl1+MN+jYiAi0ayuqYVXl+7k2blbmLtpL0kJxlmj+vDlyUWcOqyXZlyUwxLTxcDMHgYuAMqcc2PaLJ8C/B5IBB50zv13m9eeVzEQEb9ZV1bDs3O38LcF29hT20S/nDSuOWEQ0yYP0rkIckhivRicCtQAj+0vBmaWCKwBvgiUAnOBqc65Fa2vqxiIiG81tQR5a+Uunpy9mQ/XVZCRksiXJxVx40lDNIGSHJSYnhLZOfe+mQ0+YPFkYJ1zbgOAmT0DXAysOJj3NLNbgFsABg4cGLmwIiLdIKX1pMTzxhayYns1D36wgcc/3syjH23ivLGF3HLqUMYNyPU6psSBBK8DtNEf2Nrm81Kgv5nlm9n9wAQz+1FH3+ycm+6cK3HOlfTq1aurs4qIdJnR/Xrw2y+P54M7zuDmU4fy3urdXHTvh1z5wMe8t2Y30X5umMS2aLqJUntn2zjnXAVwa3eHERHxWmFOOj86dxS3nVHMs3O38vCsjVz38BwmDszlu2cN55RhBbrcUSIumkYMSoGiNp8PALZ7lEVEJGpkpyXz1VOGMvP2M/jVpWPZVd3IVx6ew+X3f8wHazWCIJEVTcVgLjDMzIaYWQpwFfCyx5lERKJGSlICVx83kHd+cBr/cckYtlfWc+1Dc7ji/o+ZtbZcBUEiwpNiYGZPAx8DI8ys1Mxucs61ALcBrwMrgeecc8u9yCciEs1SkxK55vhBzLz9dH55yRi2VdZzzUOzmfrnT1hSWul1PIlxUT/B0aEwswuBC4uLi29eu3at13FERLpFY0uAp2dv4Z531rGntomLjunH7eeM0GWOcSam5zHoaprHQETi0b6GZu5/bz0PfrAR5+C6Ewdx2xnDNFFSnIhUMYimcwxEROQIZKclc/s5I5l5++lcPL4fD87ayKm/eZeHZ22kORD0Op7ECBUDERGfKcxJ5zdXHMOr3z6FcQNy+MUrKzj/ng/4aH2519EkBqgYiIj41KjCHjx242SmX3ssdU0Brv7zbG57agE7quq9jiZRTMVARMTHzIyzj+7LW98/je+dNZw3V+zizP95jz++u47GloDX8SQKqRiIiMSBtOREvnPWMN76/mmcOryA37y+mim/0+EF+SxfFQMzu9DMpldVVXkdRUQkKhXlZfDAtSU8euNkgs5x9Z9n84O/LmZvbZPX0SRK6HJFEZE41dAc4J631zL9/Q30SE/mp+eP4tIJ/XX/hRilyxVFROSIpCUn8sMpI3nl2yczOD+D7z+3mGsems3milqvo4mHVAxEROLcyL49eP7WE/nlJWNYsrWKc373Pg/N2kgg6L8RZfl8KgYiIkJCgnHt8YN48/unceJRBfzylRVc+cDHrCur8TqadDMVAxERCeubk8ZD15Xwuy+PZ/3uGs675wP+NHM9LZo5MW6oGIiIyKeYGZdM6M8b3zuVM0f05u7XVnHZ/R+zfrdGD+KBr4qBLlcUEYmc3tlp3H/tsdx79QQ2V9Ry/j0f8MiHGwnq3ANf0+WKIiLyucqqG7jjb0t4d/VuTi4u4NeXj6NfbrrXsaQNXa4oIiLdpnePNB6+fhK/unQsC7bs5Zzfvc+LC7fhx18u452KgYiIHBQz4+rjBvLP75zCiD7ZfPfZRXzzqQXs0ayJvqJiICIih2RQfibPfu0E7pgykjdX7OKc373PO6t2eR1LIkTFQEREDlligvH104/ipW+eTH5mCjc+Mo8f/X0p9U26Y2OsUzEQEZHDNrpfD1667SS+dupQnp6zhYvuncXqnfu8jiVHQMVARESOSGpSIj86bxSP3TiZvXVNXHTvLJ6avUUnJsYoFQMREYmIU4f34tXvnMLkIXn8+IWl3PbUQqrqm72OJYfIV8VAExyJiHird3Yaj94wmTvPHcnry3dy/j0fsHDLXq9jySHwVTFwzs1wzt2Sk5PjdRQRkbiVkGDcetpRPHfrCQBc+cDHPDl7s8ep5GD5qhiIiEj0mDiwJ//41imcVFzAT15Yxo/+vpTGFl21EO1UDEREpMvkZCTz0HWT+MbpR/H0nC1Mnf4JZdUNXseSTqgYiIhIl0pMMH44ZSR/vHoiK3fs44I/zGLZNp0LFq1UDEREpFucP66QF755IsmJCVw1/RM+Xl/hdSRph4pBvGnRnOYi4p2RfXvw/NdPoDAnjev+Moe3Vmgq5WijYhBPGmvgsYtg5n97nURE4lhhTjrPfe0ERvXN5htPLdDIQZRRMYgXjTXw5BWwdQ4UDPc6jYjEuZ6ZKTx642QG5WXwtcfn8cHa3Zz9f++xZpemU/aaikE8aKiGJy6DrbPhsj/DmC95nUhEhNyMFB68roSWQJCbHp3H2rIabvjLXOqaWryOFtd8VQw082E76ivh8Uth2zy4/GEYc5nXiUREwgblZzKkIIumliDOQXlNIz98fonXseKar4qBZj48QG05PHoB7FgMVz4GR1/idSIRkU95bu5WNpTXhj9vbAny9soynpu71cNU8c1XxUDaqN4Bj5wP5Wvh6mdg5PleJxIR+Yy7X1tFffOnZ0Osbw5w92urPEokKgZ+tGcjPHwOVJXCtOeh+CyvE4mItOuOKSNJT0781LL05ETuPHekR4lExcBvylbBX86Fhir4yssw5BSvE4mIdOjKSUWcOao3qUmh3VFqUgJfGNWbK0qKPE4Wv1QM/GTb/FApcEG44VUYcKzXiUREPtdvLh9HQVYKBhRkpfLry8d5HSmuqRj4xcYP4NGLIDULbvgn9Dna60QiIgclIyWJv9wwmWF9svjLDZPISEnyOlJc09r3g1Wvwl+vh7whcO0L0KOf14lERA7J8D7ZvPG907yOIWjEIPYtfhaevQb6jgmNFKgUiIjIEVAxiGWzH4AXboHBJ8FXXoKMPK8TiYhIjNOhhFjkHLz3a5j5Kxh5AVz2ECSneZ1KRER8QMUg1gSD8MZP4JP74Jir4aI/QKL+GUVEJDK0R4klwQC8/G1Y9AQc93U451eQoKNBIiISOb7aq/j6JkotTfD8jaFScNqdMOW/VApERCTifLVn8e1NlJrr4dlpsOJFOPs/4YwfgZnXqURExId0KCHaNdXC01eFJjC64HdQcoPXiURExMdUDKJZ4z548grYOhsufQCO+bLXiURExOdUDKJVQzU8cVno/geXPQRjvuR1IhERiQO+Oscg5pWthD8eD6XzQqVg+wK44hGVAhER6TYaMYgWTbWhwwZVpfCX8yDQDFc+AqMv8jqZiIjEEY0YRIuXvgm1ZYCDQGPolsmjL/Y6lYiIxBkVg2iw4AlY8zq0NP5r2a7loeUiIiLdSMUgGrx9FzTXfXpZc11ouYiISDdSMYgGX7gLkjM+vSw5A876uSdxREQkfqkYRIOJ18DwcyCp9Q6JSWkwfApMmOZtLhERiTsqBtHi4j9CZi/AQo8X3+t1IhERiUMqBtEiJROm/RV6jQw9pmR6nUhEROKQ5jGIJr1HwTc/8TqFiIjEMY0YiIiISJivioGZXWhm06uqqryOIiIiEpN8VQycczOcc7fk5OR4HUVERCQm+aoYiIiIyJFRMRAREZEwFQMREREJUzEQERGRMBUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRMBUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRMBUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJ81UxMLMLzWx6VVWV11FERERikq+KgXNuhnPulpycHK+jiIiIxCRfFQMRERE5MioGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImEqBiIiIhKmYiAiIiJhKgYiIiISpmIgIiIiYSoGIiIiEqZiICIiImFJXgf4PGaWCdwHNAEznXNPehxJRETEtzwZMTCzh82szMyWHbB8ipmtNrN1ZnZn6+IvAc87524GLur2sCIiInHEq0MJjwBT2i4ws0Tgj8C5wGhgqpmNBgYAW1u/LNCNGUVEROKOJ8XAOfc+sOeAxZOBdc65Dc65JuAZ4GKglFA5AJ0TISIi0qWi6RyD/vxrZABCheA44B7gXjM7H5jR0Teb2S3ALa2fNh54mEIirgAo9zpEHNB67npax11P67h7jIjEm0RTMbB2ljnnXC1ww+d9s3NuOjAdwMzmOedKIpxP2tA67h5az11P67jraR13DzObF4n3iaah+VKgqM3nA4DtHmURERGJS9FUDOYCw8xsiJmlAFcBL3ucSUREJK54dbni08DHwAgzKzWzm5xzLcBtwOvASuA559zyw/wjpkcoqnRM67h7aD13Pa3jrqd13D0isp7NOReJ9xEREREfiKZDCSIiIuKxmCwG7c2caGZ5Zvamma1tfezZ5rUftc6muNrMzvEmdWzpYB3fZWbbzGxR68d5bV7TOj5EZlZkZu+a2UozW25m32ldrm05QjpZx9qWI8TM0sxsjpktbl3HP29dru04gjpZz5Hflp1zMfcBnApMBJa1WfZr4M7W53cCd7c+Hw0sBlKBIcB6INHrv0O0f3Swju8CftDO12odH946LgQmtj7PBta0rktty12/jrUtR24dG5DV+jwZmA0cr+2429ZzxLflmBwxcO3PnHgx8Gjr80eBS9osf8Y51+ic2wisIzTLonSig3XcEa3jw+Cc2+GcW9D6fB+hk277o205YjpZxx3ROj5ELqSm9dPk1g+HtuOI6mQ9d+Sw13NMFoMO9HHO7YDQDwOgd+vy9mZU7OwHg3TuNjNb0nqoYf/QoNbxETKzwcAEQr8FaFvuAgesY9C2HDFmlmhmi4Ay4E3nnLbjLtDBeoYIb8t+KgYdaXdGxW5P4Q9/Ao4CxgM7gP9tXa51fATMLAv4G/Bd51x1Z1/azjKt54PQzjrWthxBzrmAc248oYnpJpvZmE6+XOv4MHWwniO+LfupGOwys0KA1sey1uWaUTFCnHO7WjfMIPBn/jUspXV8mMwsmdAO60nn3N9bF2tbjqD21rG25a7hnKsEZhK6e6624y7Sdj13xbbsp2LwMnBd6/PrgJfaLL/KzFLNbAgwDJjjQb6Yt/8/eatLgf1XLGgdHwYzM+AhYKVz7rdtXtK2HCEdrWNty5FjZr3MLLf1eTpwFrAKbccR1dF67optOZpuonTQLDRz4ulAgZmVAv8O/DfwnJndBGwBrgBwzi03s+eAFUAL8E3nXMCT4DGkg3V8upmNJzQctQn4GmgdH4GTgGuBpa3HDQF+jLblSOpoHU/VthwxhcCjZpZI6JfN55xzr5jZx2g7jqSO1vPjkd6WNfOhiIiIhPnpUIKIiIgcIRUDERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRMBUDEemUmd3c5l7vwTbPf9vO1z5gZicdsKymzfPzzGytmQ3sjuwicug0wZGIHBQz6w985Jwb1MnXLAKObTvDmpnVOOeyzOwLwHTgbOfc+i4PLCKHJSanRBYRT4wBlnb0opmNAta0N+2qmZ1C6AYv56kUiEQ3FQMROVhj+dcNWtpzLvBaO8tTCd1A53Tn3KquCCYikaNzDETkYHU6YgCcQ/vFoBn4CLipK0KJSGSpGIjIwepwxMDMMoBc51x793sPAlcCk8zsx12YT0QiQIcSRORzmVkCofu5d3Qo4Azg3Y6+3zlXZ2YXAB+Y2S7n3ENdEFNEIqDTEQMzu771PzNm9gszSz+SP8zM8szskdb7R7d3qdPtZnaPmX3vSP4cEYm4YqDUOdfYwesdnV8Q5pzbA0wBfmpmF0c4n0hEtN3vRfA9/83M7m29nNfaLO9jZs+Z2X1mdms731doZhvMbEwH73u6md12JBna83kjBicDGa3vMRBINLNHgC3AKOD91uVJzrnvmdnXgeFAT+CnwJmEfpi8A+EfDNe3hvyrmSU454Ktn08ATgJWAzsO/AsBg4Aq59zPDmYliEjkOOfWAKM7+ZITgXYLvXMuq83zrcCQyKYTiai2+71G4HwgHfgboX3fXcAaYLJzboqZrQAeIHSo7TvAZbTZ75lZCjDROTetdSd+MvBBmz9rhnPucTN73swecs41t8lyO/DXAwO2/vI8CMgB5pvZIODfAAPWA3uBPc65GWb2l9ZcHWX4jM87x2AW8JRz7pUDlt8P/AcwxDl3O1BkZlnAV4Cq1lATnXOP7V85B/ylTgFW7S8FrUYAK51zdwDnHzA60ReYB9zzOXlFxAPOuYkH/EATiVVt93vfBioJ/bI6GbgZuAP4BZDc+vXbnXO/B/4BXNTOfi8f2N36fDMwoM1rrwITzex/Cf1Cnb//BTO7gVAZqW8n46nOue/yr1G6b7R+XQWhgvI34FIzywZagMxOMnzG540YBDtYXg30aH0M/z2Abc65uzp7QzM7HbgQ+MEBL5USGn0AqCN0idP+FXIHMAn4i5ld7ZyrRkREJPLa7vcSgP9wzrUAtB4Cd60f++3fjyYfsHy/CqCg9flAYMn+F5xz9bSOtJnZS0BZm++bDBwDHE+oMHyrzWtNrY/7D+0lAI8758LvbWYOuA74e2cZ2vN5xWAx8BMz+9yTFJ1z+8xsjpn9gVBJeBgYB2x1zr3dGrQP8BzwAvCn1uGQC1r/cjOAqa0rfqdzrrLN2/+w9S+1h1BpEBER6Qpt93v3AA+a2R5Co9Z/Bu4mNFy/f4Qs38x+RegQ2VfN7Hra7Pecc01mtsDMfk/oF977zOwKQvu9t4E/AInAo865YOs+8L+cc18HMLO7gOdbnz/unLsW+NDMfgQcBSwC7gV+ZWY7gH3OuZ8TGjW4FxjunGs5MENnK0BTIouIiBwiM3veOXf5/kev80SSioGIiIiEaYIjERERCVMxEBERkTAVAxEREQlTMRAREZEwFQMREREJUzEQERGRsP8Pk2wM3szc5mIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 504x432 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
"source": [
"import timeit\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt \n",
"import pandas\n",
"import teqp\n",
"\n",
"def get_critical_curve(ipure):\n",
" \"\"\" Return curve as pandas DataFrame \"\"\"\n",
" names = ['Nitrogen', 'Ethane']\n",
" model = teqp.build_multifluid_model(names, teqp.get_datapath())\n",
" T0 = model.get_Tcvec()[ipure]\n",
" rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2)\n",
" rho0[1-ipure] = 0\n",
" o = teqp.TCABOptions() \n",
" o.init_dt = 1.0 # step in the parameter\n",
" o.rel_err = 1e-8\n",
" o.abs_err = 1e-5\n",
" o.integration_order = 5\n",
" o.calc_stability = True\n",
" o.polish = True\n",
" curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, '', o)\n",
" df = pandas.DataFrame(curveJSON)\n",
" rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']\n",
" df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot\n",
" return df\n",
"\n",
"fig, ax = plt.subplots(1,1,figsize=(7, 6))\n",
"tic = timeit.default_timer()\n",
"for ipure in [1,0]:\n",
" df = get_critical_curve(ipure)\n",
" first_unstable = np.argmax(~df['locally stable'])\n",
" df = df.iloc[0:(first_unstable if first_unstable else len(df))]\n",
" line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')\n",
" plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd', \n",
" color=line.get_color())\n",
"\n",
"elap = timeit.default_timer()-tic\n",
"plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa',\n",
" xlim=(100, 350), ylim=(1, 1e3))\n",
"plt.yscale('log')\n",
"plt.tight_layout(pad=0.2)\n",
"plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
"plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7);"
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
},
{
"cell_type": "markdown",
"id": "6120531f",
"metadata": {},
"source": [
"And now for something a bit more interesting: ethane + alkane critical curves"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "92020440",
"metadata": {},
"outputs": [],
"source": [
"import timeit\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt \n",
"import pandas\n",
"import teqp\n",
"\n",
"def get_critical_curve(names, ipure):\n",
" \"\"\" Return curve as pandas DataFrame \"\"\"\n",
" model = teqp.build_multifluid_model(names, teqp.get_datapath())\n",
" T0 = model.get_Tcvec()[ipure]\n",
" rho0 = np.array([1.0/model.get_vcvec()[ipure]]*2)\n",
" rho0[1-ipure] = 0\n",
" o = teqp.TCABOptions() \n",
"# print(dir(o))\n",
" o.init_dt = 1.0 # step in the parameter\n",
" o.rel_err = 1e-6 # relative error on the step\n",
" o.abs_err = 1e-6 # absolute error on the step\n",
" o.max_dt = 100 # cap the size of the allowed step\n",
" o.calc_stability = True\n",
" o.polish = True\n",
" curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, '', o)\n",
" df = pandas.DataFrame(curveJSON)\n",
" rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']\n",
" df['z0 / mole frac.'] = df['rho0 / mol/m^3']/rhotot\n",
" return df\n",
"\n",
"fig, ax = plt.subplots(1,1,figsize=(7, 6))\n",
"tic = timeit.default_timer()\n",
"name0 = 'ETHANE'\n",
"for othername in ['METHANE','PROPANE','BUTANE','PENTANE','HEXANE']:\n",
" for ipure in [1]:\n",
" df = get_critical_curve([name0, othername], ipure)\n",
" line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')\n",
" plt.plot(df['T / K'].iloc[0], df['p / Pa'].iloc[0]/1e6, 'd', \n",
" color=line.get_color())\n",
"\n",
"elap = timeit.default_timer()-tic\n",
"plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa')#,xlim=(100, 350), ylim=(1, 1e3))\n",
"plt.tight_layout(pad=0.2)\n",
"plt.gcf().text(0,0,f'time: {elap:0.1f} s', ha='left', va='bottom', fontsize=7)\n",
"plt.gcf().text(1,0,f'teqp: {teqp.__version__}', ha='right', va='bottom', fontsize=7);"
]
}
],
"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.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}