diff --git a/dev/fluids/Air.json b/dev/fluids/Air.json index 63533cda..79c07a25 100644 --- a/dev/fluids/Air.json +++ b/dev/fluids/Air.json @@ -183,13 +183,16 @@ "c": [ 0.6666666666666666 ], + "d": [ + 1 + ], "n": [ -0.197938904 ], "t": [ - 87.31279 + -87.31279 ], - "type": "IdealGasHelmholtzPlanckEinstein2" + "type": "IdealGasHelmholtzPlanckEinsteinGeneralized" } ], "alphar": [ diff --git a/dev/fluids/Fluorine.json b/dev/fluids/Fluorine.json index b1a9d1cd..6eb726da 100644 --- a/dev/fluids/Fluorine.json +++ b/dev/fluids/Fluorine.json @@ -174,15 +174,18 @@ }, { "c": [ + 1 + ], + "d": [ -1 ], "n": [ 1.012767 ], "t": [ - 8.9057501 + -8.9057501 ], - "type": "IdealGasHelmholtzPlanckEinstein2" + "type": "IdealGasHelmholtzPlanckEinsteinGeneralized" }, { "a1": 0, diff --git a/dev/fluids/Oxygen.json b/dev/fluids/Oxygen.json index 74728d70..4cd9bff0 100644 --- a/dev/fluids/Oxygen.json +++ b/dev/fluids/Oxygen.json @@ -193,13 +193,16 @@ "c": [ 0.6666666666666666 ], + "d": [ + 1 + ], "n": [ 0.944365 ], "t": [ - 74.91476960299131 + -74.91476960299131 ], - "type": "IdealGasHelmholtzPlanckEinstein2" + "type": "IdealGasHelmholtzPlanckEinsteinGeneralized" } ], "alphar": [ diff --git a/dev/fluids/n-Pentane.json b/dev/fluids/n-Pentane.json index ab17afa1..879cf534 100644 --- a/dev/fluids/n-Pentane.json +++ b/dev/fluids/n-Pentane.json @@ -168,7 +168,7 @@ "type": "IdealGasHelmholtzLogTau" }, { - "T0": 309.213623675, + "T0": 298, "Tc": 469.7, "c": [ 4, @@ -180,7 +180,7 @@ "type": "IdealGasHelmholtzCP0AlyLee" }, { - "T0": 309.213623675, + "T0": 298, "Tc": 469.7, "c": [ 0, diff --git a/doc/notebooks/CP0.ipynb b/doc/notebooks/CP0.ipynb index 2395f8c6..ec7b4851 100644 --- a/doc/notebooks/CP0.ipynb +++ b/doc/notebooks/CP0.ipynb @@ -12,7 +12,28 @@ "level": 1, "metadata": {}, "source": [ - "Conversion of Aly-Lee form to Planck-Einstein" + "Derivation and Conversion of Planck-Einstein Terms" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Start with a generalized Helmholtz energy from \n", + "$$\\alpha^0 = a_k\\log\\left[c_k + d_k\\exp\\left(\\frac{b_k\\tau}{T_c}\\right)\\right]$$\n", + "\n", + "The ideal-gas specific heat can be obtained from \n", + "$$\\frac{c_p^0}{R} = -\\tau^2\\alpha^0_{\\tau\\tau}$$\n", + "\n", + "After some math shown below (and some manual grouping):\n", + "$$\\frac{c_p^0}{R} = - \\frac{a_{k} (b_{k}\\tau/T_c)^{2} c_k d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}}{ \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)^2} $$\n", + "\n", + "Two important identities:\n", + "$$ \\left(\\frac{1}{\\cosh x}\\right) = \\frac{2e^{-x}}{1+e^{-2x}}$$\n", + "$$ \\left(\\frac{1}{\\sinh x}\\right) = \\frac{2e^{-x}}{1-e^{-2x}}$$\n", + "which yield\n", + "$$ \\left(\\frac{x}{\\cosh x}\\right)^2 = \\frac{4x^2e^{-2x}}{(1+e^{-2x})^2} = \\frac{(-2x)^2e^{-2x}}{(1+e^{-2x})^2}$$\n", + "$$ \\left(\\frac{x}{\\sinh x}\\right)^2 = \\frac{4x^2e^{-2x}}{(1-e^{-2x})^2} = \\frac{(-2x)^2e^{-2x}}{(1-e^{-2x})^2}$$" ] }, { @@ -35,58 +56,155 @@ ] } ], - "prompt_number": 2 + "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ + "a_k,b_k,c_k,d_k,tau,T_c = symbols('a_k,b_k,c_k,d_k,tau,T_c',real=True)\n", "alpha0 = a_k*log(c_k + d_k*exp(b_k*tau/T_c))\n", - "display(diff(alpha0, tau, 2)*(-tau**2).together().ratsimp())" + "print 'No derivatives'\n", + "display(alpha0)\n", + "print 'First partial w.r.t. tau'\n", + "display(diff(alpha0, tau, 1).simplify())\n", + "print 'Second partial w.r.t. tau'\n", + "display(diff(alpha0, tau, 2).simplify())\n", + "print latex(diff(alpha0, tau, 2).simplify()*(-tau**2))\n", + "print 'Third partial w.r.t. tau'\n", + "display(diff(alpha0, tau, 3).simplify())" ], "language": "python", "metadata": {}, "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "No derivatives\n" + ] + }, { "latex": [ - "$$- \\frac{a_{k} b_{k}^{2} d_{k} \\tau^{2} e^{\\frac{b_{k} \\tau}{T_{c}}}}{T_{c}^{2} \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)} \\left(- \\frac{d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}}{c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}} + 1\\right)$$" + "$$a_{k} \\log{\\left (c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}} \\right )}$$" ], "metadata": {}, "output_type": "display_data", - "png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABMBAMAAACPE1EUAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMiKriXZE77tm\nmVQzv5s6AAAJ0klEQVRoBc1ba4gkVxU+VdPv6ReKjySu3faaB7JxJ+al/rGJBGRRpjEkQv5MIxJH\nMWw7KruIyzYxrBuzy7aJIEZlRpOwSVzdSX6EIK7TRpQlRLeNiQYUpkB8EcjM+txNdjOe+6p7btWt\nqumanuj9Ufec75zz3dt1q29V9TcDsNVWfGirDLx+QjTjzWVpvPSo7AnRRNHb8GzTho6NTYhmrHGr\nu57n+Zl3PXiNrbDct6EhLImGFnyHOluwK8OPQq4JUC0OrCx/tKIhMImGFkz3qZfedjsvQXWE9fm6\njSS7yW9xAo1B7TxquKkdFz4FFQ/LXSvFVN8Kh8AEGjP/gY7pR3kHTllOZ/73v5b5xf0eTF3yS4DL\nrAQLVjQMJtCYBZWm6Ud4+Q67koPtC/BxApV+/jbiUTN7nnoJdjRNoNB5LQDY3alO5l/hyPUwSxbt\nvkPhDIFML0dFLHg0TTD5i0HA6pcGtjP3KMwNdPoR9+3aMay1tuHGO9E0wbpGP4jY/epZG75GvgEv\nlO+xpSDWisCtcDRNMD1vnRLNuuxB5olNovijEQ3Be6W3QZqRwJzMhRCkgW/9QNtAWDYIHGFm43hF\nzRWsOyDsmid6cczrD1K+8uDvrqQxbefWtR22zC9bDE24NHk1GXlZTnKxThlwN1StAnvgrcox+0bX\n9A2vaK56DI1Rxp1jnTBmIJz8G/B1Dn6Ohoq9Ylv5VfgFXKIcsz8zMH3Dy88YbgyNkcedRi+MGQgj\nz7z75cc5eDUN3fTyS3Xtz2vTtFZM1/TckelDJE0gD113KYwhkv3k3deIgDt/1CttbKxzTF6ZBfxC\nnYOTG/hlcvYf7fDE66w8CP49KgBw2/5ZT0STacIs+X+GMUQe6EyvQ+E4Wo13gLwdIXZWYPtu/sSt\nYsJwfyfbYwwOfrbql0bMNFvxP6ZPvOLDsFKHr7L6RBpSpszqOWXR3jkBlWWA9yGGe7cYnGFNjmU8\nwCBv2R/Cjjazik2APwObRaDlox8DLh3Cc7jvDnBBE2kCrMx1XrWAUDnLdnKHff2fA0dMnWE9gUGm\nK6tyJw7/lJsLNwK8U4JG55p7CI2t1gGrGnWAZBpap+zHlEH70jI02lCcQexqKL7GF5VhA4FBfiiz\na11psO5y6BBPmjVGYm/fBweX5H4MJtPYGPZahgPcd9agnh99qAsX+BIM8OxoDKbakqrWI5wfOT0k\nnjTnlsOYQHA1q+sOHMgcgmQaG8nBoQVtePBEteP+6mtn4TzgJckWlWAwp2pKOPW+cqz92sgKI4iP\nrdPNP8GuvzXqyTQ2kjM4eKi5/fINeaixl+XPwCmxqASDL6uKQhN2dJRj7Y/1rTADr4KfdXtwsQ2Q\nTGMjWexaUGf+87fdA3vu9QB2HBriIwwuKsHIjePoEQzHtdXo+Ju/cvP8IPPIT7A8kcY2xFzThnLs\n3sxSWwTZogpLY8JPPK4M4lMKy3sy8RnhqCsmU1sKhyTyVGamJ0y2qKJpTCEJ/e5OfMJ0/4P5+IxQ\nNLsip74eCingcphvc1ssKjd9TCUl9SfFMJFpbid3d2TQDLxJuM6pvYKzFH3L8AvTLKoq3lDG1vu/\nKopjYupTlhdmlaL6FIuqSuF/PPUxFtWfsjK2cequ/dFRjbzV3tnOqeNjxDY255XJkQev9enz8IZr\nWcPnPvpKPgGbTVpO3dnJhnhPO90Q8tNbpj658xJmmtgFs9Bqvb/VEm+YcofBs76t7f//Wo+UsNJN\n3UoXvGDid5jycLNLskQSna520k0dKJ0iC059KnZzvEuVJfXZJs24TztJDwI6k1omnYyoqf/25A1D\nBsXeTYvflVWJHZOwck2Amx67Fu/Sha5f8AFx0/b9zRmKzhDW1NQVRWldWZb+0rYFtEJMwmJykvdN\nHt7lJ620fXMMQ9EZwtq3AwS1uKnfGUiOdpmEVfEw/lmeM+ef61UGjt18OruwJvgaM9G8mU08Vcpq\nJmFxOanJgVxX4hDzgqdSLL1P51qCCppdVla4D/0iGE5RCJOwmA7kCDr9odd6KmWcXtFFCGuCaq0b\nplTa3FwnHItBqA50lcqb7Spr3J7SWWttC6q0ud9YKyJBqgM9o7IaTWWN21M6a+1qOwwrbe6RcCgO\noTqQr5BNrceVxMUonTVvN0Gzrzy599U7dvekNpeNvVsFhSXCgyb7UYi3SiyJTGKdISsRPNr8Bwnl\ne1CagcqQ/76LMkz8Cat5pJQ+IzOY3Z54q16QRmInf7zneUE6W3HmIkFzAHMjKA9AaHOVGRKTptPx\nscW6bzLj6X0XF3o+UliXZvbfPmYxCF1AVkpWxAr0s34Y4IwHhToIbW5qOTxYte9jhrAEMMosOTpY\n9d/WuXjpFwUMQseUH9KSFTFDGhwA7K1DEUBoc7URoZImGcsQlvARoNx1kEC2jH9KDnYUZukJXeAm\nkqyIlZomIT9HUpur9Xls3+G2zlFjodikJveWwy/weK6v08DxL/G5HoHRtNLZZKV5sy7kLfYNKMsv\nfanNNYYsdovnNrmqxBPV1LWwVLwePsZD/o/u3DvHj3hwZ5TFezudRVaC64y6sIMXCG1lfiq5Nodf\n2SELHYecx1UlnienroUluN1znuehBh7/Ii41tISKg0bZ3B2tdOw3/KCsFCms8cHwQDcYdHP0gWvW\nQwQlJHyv92E5dSIsffr0EfHx8f2i2vV/+tTE/odAogg6WA3LSsVmhLDGaLAFTgnb1nWbHaJdWsZD\nUcD5Vmvn91qtEUeVsKQ2v/zOHrhtTBZNT/1gXWHY2+nAIitFCmuSzW0SWjRxW9eNXzBczuaqEg/I\ns06EJXpS3U5WlWvYv7GykJ0uSlayCmtyhMW2NGSH27pu/I9lakPc97iqxANq6lpsOo7vc6om84cj\nyvR3GCj7VxvG7HRRspJVWJMjBPZmph7jn/wd7vNwrYddpYuiLleVOCanToSlu8ARmyOPq0NWbZ0I\n3KhA7O10KWSlwB8X5Z/YuKKPgnRmiQ9W6rJu37N9oSpxTE6dCEuFF5/lEfOQOa/92Y627XQwvqxU\nGhJSZVaWnTa3p5sK0gqSnLoIaFgl+j09KeUZHw7WGXRjKRB3BDi5W/Ikml+SBmgFqThUGPYaJqAw\ndS36dwbCus6gG0eBqAZPBx+iptZXn7kIBSkCZizTXc4lDq5HHDQj6sZRIG5XkzSY3QGMOKCfQ4z4\nphy2k/hNXoC+PwHjVitH9sXTdRF40hrfFGh8NTdVMdmkZ9LTLaQvnUgl/kVD2nYibeGE6grNtETV\npbSVk6p7KC3R1DBt5aTqnk5LlLow7YChurT/EZF9OET1egNO0ktWxITe6EUEXkf4FutNK3ECP07M\n2P6EbDvNGORHjTTlm6/5L9dwOhTAi9ZrAAAAAElFTkSuQmCC\n", + "png": "iVBORw0KGgoAAAANSUhEUgAAALYAAAAoBAMAAAChsh0EAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIquJdjLdEETvu2aZ\nVM0GsGrEAAADx0lEQVRIDa1WXYgbVRg9+ZlMNr+DSosvZqBFEKQGd6tvy0C3IC2WEUToixsVrb6l\nFAoFIQERLYXullZE1ocgfRFR8tAX6Q8Rf0AWYUARH0RTLPqgLLVr2UK3rN93586dO9OZTHbpB7lz\nvnPOd/Zm7iRZ4EHVafdBJd2fM3M/lcVUh1kOqT+m+0xb79LwkTQhzv92xb8pL3/0Rhc4G5cT+vKH\nCWQitZobYx8p3vssF21eJ9epqU/oen2AXznsWZH49+RcVp/MtkjHON/HD4zHgmjRjZlcxsZkXVPf\neR44dmYEcyDIvK1pibAxSqTTyJXdw0DK3lXLCrzTXMs37NB2I4TJ6LtkOoUtDv/ylPSlQingQgqf\nTNesS32ltELInLGvrSQBcrejfbQrfzoKieZWWDeZnollnbRDM6OqcEU5rWt6WgNcm92cDwPz44iI\nih3ta51oz93xkFrqhpjQyFg2HcUUY/uKZ5cGyqrATwrhmRAy8qq22VdU/V8FBYhnN0dRnTst+5+Y\nmnc0wrgH5PYcPBZQlG3OHeijvGd1YUBk0xHK7IIrrmIJsmmOpkWdWPhKXEuaDeZd4LRVuYnieSFS\n9ksWLmKpWxo6xPSGTO/yGmNpoC7IprkNf678Ch5mH3r0+vGKgLT8B3MFtQHwpmAo+zVg0X0c+Q4T\nLZF9HnlPGoiT2Tw39uce9cxv2c1f2nW7ICAt66htoOTA9L+TKjbfpWb7e5m96JGHHAgMhGQ2z7V9\n+olLh7skoPBnGw2XkahNzAzQc1HuiLZi1+9Q9qDlnuJULA5pIQcCw/za2ltra3yGPEcH0yF4i15B\nNaxcADfRa+N3dAujF23i5L5Hj8ztFw5xT8hBe/INhOS+tbl14fUXY/Vw0K2j5+GLutW4+h6/84qN\nj4Gl/i9S7zkEmrR5TxqoDbLDOXoOinJAv9xFw6m+WkDTP4zKAEdd/IGz774gXE3ecs3G2640UCuz\ntbkjMP1nUE9G7h7Mvfsfeg5HD3kkGE/dsnJfP91HY2uL9k831eZ19rIjDdzJbG2u+M1lFmJl0MmJ\nOmQsu5r2OY7/7FBfGQdkaJDZvhDSgVFd6/6jB1w0Ovz2ZdEvNqr0QmFZMprhTEDxNTqnK+HsdezV\n9m2+Tlv2yKn+NqIGFZJCs07PRWKdmDsgDpe/E3Za/HhNrE8mqhPFRWuiDGT+oKbPz6dLvkKfrJ3W\nStZgcZzlSNPry2mK4j9QaJuglHWU9Nu9zUhln2Kw4ij3tkDuXLadP0Y7qZPTPAS7rJ1E47NppnLu\nNK64R/yf8j+av+zHorG6yAAAAABJRU5ErkJggg==\n", "text": [ - " \u239b b_k\u22c5\u03c4 \u239e \n", - " \u239c \u2500\u2500\u2500\u2500\u2500 \u239f b_k\u22c5\u03c4 \n", - " \u239c T_c \u239f \u2500\u2500\u2500\u2500\u2500 \n", - " 2 2 \u239c d_k\u22c5\u212f \u239f T_c \n", - "-a_k\u22c5b_k \u22c5d_k\u22c5\u03c4 \u22c5\u239c- \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + 1\u239f\u22c5\u212f \n", - " \u239c b_k\u22c5\u03c4 \u239f \n", - " \u239c \u2500\u2500\u2500\u2500\u2500 \u239f \n", - " \u239c T_c \u239f \n", - " \u239d c_k + d_k\u22c5\u212f \u23a0 \n", - "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", - " \u239b b_k\u22c5\u03c4\u239e \n", - " \u239c \u2500\u2500\u2500\u2500\u2500\u239f \n", - " 2 \u239c T_c \u239f \n", - " T_c \u22c5\u239dc_k + d_k\u22c5\u212f \u23a0 " + " \u239b b_k\u22c5\u03c4\u239e\n", + " \u239c \u2500\u2500\u2500\u2500\u2500\u239f\n", + " \u239c T_c \u239f\n", + "a_k\u22c5log\u239dc_k + d_k\u22c5\u212f \u23a0" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "First partial w.r.t. tau\n" + ] + }, + { + "latex": [ + "$$\\frac{a_{k} b_{k} d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}}{T_{c} \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)}$$" + ], + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAJsAAABMBAMAAACIdNgoAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIquJdjLdEETvu2aZ\nVM0GsGrEAAAFE0lEQVRYCe1YXWgcVRT+Zncns5vd2QwFfTSL5kHwoQuJiChxiSnUJsgKIgTFrBZ/\n8GkQKbVGdp+q9aFsrRY0KkH6pqarFEVa2xV/XoS66ov4ku2PFUmJKRFSbUI9587O7J3snZ+tPvTB\nA7nzne989+yd2bsz3wToO9Jv9j0ldEI9tNpvUWv3OyNUby5/I+qPv/W0HSqMV8w1bkGiDbRejaeP\nUOWtX2A2SXRXhDBeOY/bkWuRth1PH6FKj7UwsPckjIUIYfxy5uOX4oujlQcnojV9KKbyL/ehjpR+\nmd0ZqfEJ9h2R0tfekZJrUkh0BDwn1zflBKdGN8aLPiY6kTuk1336pl43Sj4mMvF1SFZ8+la2YNR8\nTGSSrEiSfFNKCCZKNMQI7bYdzziy/Mg0/5A4to3NdqAxNm0xM1DmMToOWIOrSB0mYfU3XMYDRULp\nQzhmO3C/pTGDKv39dJJRaBjzyC0Az5JoycYaBmqEXmzgDC2IoPYu5sS6DgJmIUm18MitY6AEg7/H\nMzDWULUJfW9jNy2IYGJ+8jNukLxYRF705Sw4MguolpGukOIS0pvYz9K3YVyBgEMFzp3IW5oLA4/V\nIpZgJ5sPF3AVtNR79QlepLlqCDhU7M7Uv5vqJgGo2sKiaeW/eGUdV/iaLf9ctY1NDLZ/EDBD7UoB\nU1V0vpR9MokhfrTcgQ+BDb5AF/BpoShgqo05SzUvgDNG7ty2EzO7WsDcRAP6Gx+R8Pn7HxqpOXB6\nqhEwM4TepdfLopxamNE7OgmGzFSWjuqVoigMlh5xN5cElXNCyLMYcVaXtxI7OjoJhsz8v3QDXQHp\nsfTv4eoNdGL/6VJiOu4TcT+0LgmNgpTIkG5C8cLvuOnpoYwDtpLuJdlxJ9qA/vsRepamCr0KYozY\n7yDsuNkim+kad1rmoScSzR4qgGDHLSxyUpzPsPKsZpWsqiM7brbIyItqoqASXVKRSo4dt7DI+0RZ\n95ssZ4p2WTk1iJQt8gWFKFdRkMGUbJE/V8gyJQUZTMkWebjGOu2v97f/fesfRWfOcNk5qkefRd4q\n4Wc62ZciMhXkGk719FaRP5cNrnw7YxVvaRqB4SayNcbAt85BHg3Ly3wGF8iev++r814xtcrwUeB0\nCynboZ/zqh4wSx5MVjzIIIcZvOAxptgUNWC7jXSH3fCqHpDabbHIJj7BXk+nu1dCeiH5U1RHJ8ue\nCl67XouMka4MxlUn0bpLcqibW/m28Mei7rZTWGQ8JbWDsxRk3VXSrhGfcBiJlvDHvnZk985sscgG\nXX3znman55pzTHR/beL8yXXSjcoj3dX1WmSk28CPKHbadU6Stp0bOq+ODDK5/woNtCtXVi6+t7LS\nJKiwyOO7gV9ZJqLTjrYdzZss0ShOlgwy5eyPRXRWF2CRcRaWo6P3BhG07ehVQq/TaPDlHGrQPzWE\nPxZlt12ARX7sOMlFdL7ZYzbtyAWjzBx/Qq5ALyPCHwuVe+2iLLImvtHk4rVzJWRaYipe58PoiZLj\njwXntouyyDq9fLgxZDlo0SU8f9zdxlwKscimtxfofl9Dk+V388Dh+WOkGw4jxhCLnKx3ddrXx23O\nlsRIwPPHXQ2jEIs8WPBLOavKS+kthzG8IbaG6iO2agLyWau3kK33cjGZcYWO3ueuN+ZVEx9UkXE4\nU3lemXKcuQrNQENBIttWsTG4U2rNUTUdxWqH1Io913e2e1rqdtoTaj6C/SCoflNQIYw3ak71H3/h\nm+3fsAlYAAAAAElFTkSuQmCC\n", + "text": [ + " b_k\u22c5\u03c4 \n", + " \u2500\u2500\u2500\u2500\u2500 \n", + " T_c \n", + " a_k\u22c5b_k\u22c5d_k\u22c5\u212f \n", + "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", + " \u239b b_k\u22c5\u03c4\u239e\n", + " \u239c \u2500\u2500\u2500\u2500\u2500\u239f\n", + " \u239c T_c \u239f\n", + "T_c\u22c5\u239dc_k + d_k\u22c5\u212f \u23a0" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Second partial w.r.t. tau\n" + ] + }, + { + "latex": [ + "$$\\frac{a_{k} b_{k}^{2} c_{k} d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}}{T_{c}^{2} \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)^{2}}$$" + ], + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAKgAAABQBAMAAAB/pCPzAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIquJdjLdEETvu2aZ\nVM0GsGrEAAAFkElEQVRYCc1YXYgbVRQ+yWx2JptkMhQUfHGXug8FwUa6PoiLhnULtV00ggiL4qb+\noCjCUETaurIBsVIfZLVaxPqzlIIvdRul6EMrRmwVfKjRJy2UTVt/kNZ1RWGr7RLPuXfuzNyZzJ2Z\n6EPPw733fOecL/fMT3K/APRtxht9lyoK5xWxfkOZTr+VirrShZMsmvt1/+OKtHShYvN6yHYASkYj\nXaEq27ROQ6mFGZqtSksXM2EDFNtYY6aqGz+s2oMx1obBnccBdqch1Sy6ZErLf/icMh4ODlq5P8Ko\nhLwyIbkJnHwjcykmbZv5fExGOFxaDWMS8nlhi+Srnd37KT5YlbNeftvnd33mg1XLcxQcD2asSUDh\n/B1fnJeQGIfKC61AkiFfjiJMwzOBFJXLyl+AF+UcrS75JfgIdkqI2tHqALnfzrwvZ5kt2YfRgN/L\nzdyw2fnCMUen2vludwWzNk3WeO66sRl6K9H0sSmLLR5jo3rYaw2twMA+TJr7BZwH/9q22WGQ8Soc\nseGuCkb3WBmaQF/BL6nbWrSMNP0AFBcAnsCEJRv+5Hn7INtm0K4mnMKnrAGQeQdmaxQ1OgDfAeMn\nt6cVV+nJ1OkWnwKdkyKGO6LhGxsewRZsgOyByU8Ywe0I/MRW0UN+AeZqYNQx4yIYa6xVxHBHBL0F\nOr61e3BVHsFB2FmwxLLnPFeBJbC11n0jcBlo2w3cGTVHEO68tKLDeG4Cyv6GHzza7EkmwLk2LJYs\n89OXVuES4CWkVstY0iZIX4Ohzrdw4fs5O4+kVVEUN5vVwsMalE9i3o1wmLdaHMGbwqAf4eORClyp\nAQx0YNaKIxNxffSWdVtgemsbYHaiia8+tgqbjlU5tOPOe0cbudc/wOypbRhNZVtz8zVeQK2ylQcN\nLEznUrE5yYdy9QpfUqvMPGioer/mgKmmszDKyXirVOtCYFrZzanYgsn9thrkkfx+W5VIgs5/bjVI\neBX6vp/E/215FbbZY0vqQ12PggRQgkNdApZASoJDXaAigZvgUJeAJZQSe6gLVUhAb3UWPNQ5Jcek\nUoUz74/tcJzQoY7h+GORzGR1RmcAtNChjqGw1+Zz7EjqLNvBMxdXZ/jTgxY61DFUT6xjSZ2RkHLU\nmdnE+vChjpFmMS2ZkTpjQoqrswwdY5xDXZBgJmn3KMo2cCHlqLP1QSrPv+gtY1akzpiQctRZ9HYy\ncdJJ/iS/kDKrcgyPS19/xaBiPRhR+n4hVZCeWyq7GW5i1fkqm5IOfiEVlnLbYcYipuFaUr5wHh3B\nmWX+Prjxn/W/V96E4QYBnzE0alCrM94r1moVyNeh2GQncvS/jOLjuFKd3SNqs9hxCwoNgEcZ9KQI\neLPOrgvz1epsSSQ+gB23YcDG0zIru+KRiVWpKlag1d0lLQLqTCgiaABstMEAOM7T/+KTf/SRqtXZ\nHF5EYUwwGxWjhoB+maGuICPPJY1TZ+UqK6Yhwzp+6MxpmxxG6goyliRIY9WZT2AU2B19r9slhhxz\nXEEmke6KU2flFsunIUtfWcJytFNPkDFY7DRWnflI8TH1jLXvCTJ8jpeXf353ebmFKbHqzNc+Pqae\n6XT3XUHGcWen8eqMNJdj+Jj6jBQpBZkg47ggjVVnwx7TERtrtckqZ3gNJ0+QMUxc01h1NtPgHKAt\nds9VUY7n5jmwSJMryBgmSGPV2RJtz7Pigl7j3q0C9ASZIGURlWRxazlHvi243E/zBJnRFEGcVZLl\noC8Pl2VL+O776wkyEWKzSrLwV9xNNxvQ4s7QiAumXQT/bsycOGpzjvCvV2JurROViv8V9Gt5/7WX\nSe6W3RTeD9G5+Vp0TB15Kjpc6ETHlBFl4SFlaXTwOvex7JHzbH/960/34HKhzHZ3mWYxFH3vieaa\nNFxuru+P1H8Bv/Dea92wCuIAAAAASUVORK5CYII=\n", + "text": [ + " b_k\u22c5\u03c4 \n", + " \u2500\u2500\u2500\u2500\u2500 \n", + " 2 T_c \n", + "a_k\u22c5b_k \u22c5c_k\u22c5d_k\u22c5\u212f \n", + "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", + " 2\n", + " \u239b b_k\u22c5\u03c4\u239e \n", + " \u239c \u2500\u2500\u2500\u2500\u2500\u239f \n", + " 2 \u239c T_c \u239f \n", + "T_c \u22c5\u239dc_k + d_k\u22c5\u212f \u23a0 " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "- \\frac{a_{k} b_{k}^{2} c_{k} d_{k} \\tau^{2} e^{\\frac{b_{k} \\tau}{T_{c}}}}{T_{c}^{2} \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)^{2}}\n", + "Third partial w.r.t. tau\n" + ] + }, + { + "latex": [ + "$$\\frac{a_{k} b_{k}^{3} c_{k} d_{k} \\left(c_{k} - d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right) e^{\\frac{b_{k} \\tau}{T_{c}}}}{T_{c}^{3} \\left(c_{k}^{3} + 3 c_{k}^{2} d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}} + 3 c_{k} d_{k}^{2} e^{\\frac{2 b_{k}}{T_{c}} \\tau} + d_{k}^{3} e^{\\frac{3 b_{k}}{T_{c}} \\tau}\\right)}$$" + ], + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAABUBAMAAABTgQyuAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIquJdjLdEETvu2aZ\nVM0GsGrEAAALOklEQVR4AeVafYxcVRU/s7Ozb2bnU4gaEmNH3KQaA4x0xbQ0ddKPBGhTxtRoGgy7\nVUFFsSMSQrE6kxirVVMWK+WjqBvSRIzYjqaRkCodhGJSSBmVP5SE7AC1xrSuWyUu2tb13HPvufe+\nN+/Nx9vZbUlvsvPOPb/f79x75915M/veD+A8tvj9wYO3w4JV5xGZaDN2O6yNrN/QoS4LRpo28Q67\nA27MBS1mZ7DZ5Wjpk0eIGfvbns8CxF0yxm5+6NZyl+UWgraz28FTtXfDAK4gHa+KeawWL9wYa3yb\nM+fj6LT5RLvnk8m9DOk65qK0+kzNgjX2YSu56OFAvdshM/A+SDWQnSFFZNYSaqxpJRciHL26zTYa\na4O5JxMfbcDQ3b8GuEfmL7dgxpxJK7kAYaqWajPCqZ5GTPziq5rvfRdsTJP6HKQKyZnAkpHTgZAf\nsMv6xGeKboaNuZF+9lLNwGqp8UDID1if+ZpOJyd0SIGNuZF+9raVAqslioGQH/B08jqdjrypQwps\nzI30p3fPHqxz53FPsfhP65xZErxMpHznEebhcc5qlN5tQAuaM9k+R6+Jei2f/WyDhznMgf/xnCud\nfP0jz7xuMleaEMCD2VC/YppMzD0lgHvLXP85DnyPcfs7BN8T2Ax3GuKNJmzBbKhPsZjMtlLkjKfc\nNbp/m478gui4K5uGX8LdJjOVMzF4MAvpVxgdBxjKJR/w1DNfKmc9iLubqbv7MGL3x/ROpawLs3nz\niiPvXYe/ZEXLjGxoxEavFm9gdHQdpQBRs+feoJzGJEO9XjLKs3VGN8hzcItNqNREzx+zefOKd+aG\nZ2BQXGYqfwX1ZRi5Habh+gLmEJ2VKE6Edp/B7GHj98GBspTsyEWEEpwZ/Jl8bV2xskUR+GOKMu+D\ns5euXJ/DQlNl+JeslyjAFTBUxfkg2gQQKID8JGmMcvyyrQbHcIdWkfVD2F4SafFPyx+AloW9rAgC\nMET60lKzMFQER1yDjoGj1vJiGbsVfBFoQaIA8gqnMVzpiVuxfamEZMx+WkoG9q59gma2ChN/oUi8\nZOv4EoAJvB8tMQmVEsTHsdYpiJ+T2+RhUXkH/gm0KlFcC+0xjQmObg+Dg9/sQpLN6yTAq5BTPVpL\nAGYJ5hVWCjAF5Wj9Y3k4Q6ehiu83nZ6VsdVgoWqPGcweFrPpGQeEhDYTY588WFMhpQMwZs/3WGnA\n/nQu85tvzcKbgLte7CxcFX4FnPxTpWyhuES6jmnMHtk5B8PN35MEP09QtDEZZ8WiArBWdrhMppj8\nVBSyR1D9AXhcbhP4LThH4GwJL9MGRQKdLo25xjsOv8oXSIK3N7bnXBh1ljTwEIC1ssNlnJEPXXId\nbL4Bh9q+uga0TSB67UGI/eDneCoMiuW/L4bQmGu8O9Z8dKRKEtiwHsu0tLGqSPljLeR5JW6ITZRk\nAbGzKBqc3BxTNRndz2NYGKfEMSBNlClZ1aYvULwvNl6QpcXOojZc/HhUhYwuV32wME6JY0CaKFpr\nCxYkfhVG5BrkNhFjZHID69RYjOr31sLs6QSkifKoTVycuN02AfmbKtxE6LspnDSsqt02geF82LIQ\n6+22R+hxbGG7bQLeGxC2sEMcbXYgLDaM34hhW6IWVrlQuo2hC/85tHKhhAm+YPc8wBd6Viy0INkM\nOUJoYcjxupHt64bkw7ks55M8z6m7wm0y5/bzPG+/4SNb/LIdc8MX3FVMTPntHeftR7Duk/nBF1nO\nvgf9Fo8vsjN3MS43vINl04pcL+KbVpd6oYc6FxOhVCiaxDtvPYgfSze7ph8KN6fwDpamM9uLuBwt\nuOlOPmjGeE8nVJuHgyU9zuJuRo6tAA99V5CsazuLp8A8HCxfz7HYU9O/G2966IN5f2L3dhaPPryD\nxVmKN0GUNcZT1K9bhNNe+kk/Ht5mr/vnO2bDO1i2v7OgxR2HAWciNuGlLyn76sb8075cVzK8g2Xr\n3DiLXRUDOi/hNZmdNIoykPflnlLZ6JqlvrhIrnw8YME9OFiiR3/nqd9W3Dodmx4Tj4Z0Y662szwP\nRzXoCaK5gaYnpbo9OFg+CFd6arQVt07HRT9u12KutrM8A4Gug6Fc0J2sHhwsW2AsZ08AgMSRrXNz\nVXeeeq3TcdGftCXMTRR19piOPEGi6nWuMKEHB8uDsKTKMnkkceoTSwP2r3c6LjrVivzn0av+e/k/\nCuJJpGjGztJycY4sW9GQo0Ka9qflgsG8/R+DovHh0mX4aMNupOQ7z1lLCWWYtIkyFnQ9HV+6eAAF\n0QIkxiFVY67eWM7GoqfouyB9RqWGJGZcMB6quxupw1TNnRLKz7hTqhcXs/K2bKN1OpKj6PQBHsBT\nUYdklbnPmTL82XdyMndjA/6t0JXyeG/ZsDF6atnZVa6ZKOXgOciOu5jCPxOt65Rt40nVOP1lDqTd\nhqfjdv0o+uAMkm8CONyAQZyU5N5mCow1ZJwuyuPOqvM/GSXr8mhcMNSvxyacokTkq1LGH4RK084D\noBJ9ltxsi89QibMvcUB04Om4XT+Knha31asAV5UhjpHksp0Ft+eSmizGa8FnrWqPfQO+SRB/F0ke\nNJJ5p6piOhild4/ho/VCXM/atvEM6QLWWk5Z0xGPfI0jSNHZOLUH1Zor7SzL1pa2wJNlWdbMSDhz\n0OkS+/srP8Gn38YF85W1TxN1oEgH/WKUYtVYVCKkvPmVl7E+K0e0xgS8FqLL6QTSpXUFInQq1NRl\n7h2NTHPTmveosnpGm76IK0AXTGJubsZ2wcS3wKXENftDSlkZu6YOIIoO7hYA+Wd+PId+Pa28RQpc\nr7wWotN02tDlOUjSQwY1dWln2Q0DDVOWZ4Tn5H7x9P0KgiwXzGUNfDouWgX//mh9DoxyeQmoqDDN\nKKUQsNJt8REINrWW7uj0KB4G7N8ytO/Q4GI3MyN4qPwibgxqlgvm/QfXy+wu3Mt5fhiLLKNMPEBF\nyVIjlVSElfGmbfEhSK+lO7r8nCfGlVYcyM6SmMRIbgeITk+f+NH0dB0zbwM40CCny/UFlwvmnwiK\nFj1RgExJxtjTykgOhs9RUfzHnpVUn5Vui4+osGp6+vPT0+JjhsJKVU4nmA5yLfj1YhrtMTS4YBPb\ngRq/u3NlOPAEnUthm7JcMPL8SnImF5GBeFXK7CwMn6GiZKmRSqpvK43FRxVQe6w7uqx0uGEGl3YW\nYUtp0HYghNeCp2RrDq/KafLDVIxHZjeeRa4Re349h3ot+M4OnaaiZKmRSqpvK43FRxXgtdBAnejy\n2+JAGbXRtUVZQawvlUenGm0HyvFaDkH6DemC2YF5ywWzERx5TZYlzKtSJmswVaSiZKmRSqofqBQ1\n1Fq6okfEBSy6f+61Im6o2AR2sJGdZdmhIv7EQIcVNV5L/IWjVXK6kB/GcsEMPntIMr2vrPzeC4/h\n1wsWJUuNVFL9QKWoxNcxstt0oMfQmKZaatIpyXA/p2g7UIdnxAiQbUr22AVjMDtqUVp0U99W2LFa\ni0x1oMsf7kRNNLjIcg5oO1AnXuMcH8+qhQOwC4YR97FFadFNfbfG9L5rQnQCym8wO2XH0Qndy+Y4\n5H8q5HbgrPto/DDoL1QeGTcjqGfRabsF8VrzHejDeS3JVKEuOxU+CdZ20DQVtPXDeMlB/Tb1/SQd\n6OIqqVrk2YNlGeoFWtuBWXxs64dhUqdjm/p+0g50770DKpGcUJWs7eCt3dYP4yUH9dvU95N0oK/y\n08zDzuJXbrFye30H2uibvcCTad5N7nmGt7O46yxqb8h89O1xk0279xaJnwqY576A/AWcjtwXMLmQ\ndpaAaouSvqsRMExIO0tAtUVJ/yxwlHB2lsByCw+oG1n/Bx1oNENKOewqAAAAAElFTkSuQmCC\n", + "text": [ + " \u239b b_k\u22c5\u03c4\u239e b_k\u22c5\u03c4 \n", + " \u239c \u2500\u2500\u2500\u2500\u2500\u239f \u2500\u2500\u2500\u2500\u2500 \n", + " 3 \u239c T_c \u239f T_c \n", + " a_k\u22c5b_k \u22c5c_k\u22c5d_k\u22c5\u239dc_k - d_k\u22c5\u212f \u23a0\u22c5\u212f \n", + "\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n", + " \u239b b_k\u22c5\u03c4 2\u22c5b_k\u22c5\u03c4 3\u22c5b_k\u22c5\u03c4\u239e\n", + " \u239c \u2500\u2500\u2500\u2500\u2500 \u2500\u2500\u2500\u2500\u2500\u2500\u2500 \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u239f\n", + " 3 \u239c 3 2 T_c 2 T_c 3 T_c \u239f\n", + "T_c \u22c5\u239dc_k + 3\u22c5c_k \u22c5d_k\u22c5\u212f + 3\u22c5c_k\u22c5d_k \u22c5\u212f + d_k \u22c5\u212f \u23a0" ] } ], - "prompt_number": 4 + "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ - "And then with some manual simplification:\n", - "$$\\frac{c_p^0}{R} = - \\frac{a_{k} (b_{k}\\tau/T_c)^{2} c_k d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}}{ \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)^2} $$\n", + "Aly-Lee starts as\n", + "$$\\frac{c_p^0}{R_u} = A + B\\left(\\frac{C/T}{\\sinh(C/T)}\\right)^2 + D\\left(\\frac{E/T}{\\cosh(E/T)}\\right)^2$$\n", + "and generalized Plank-Einstein term is given by\n", + "$$\\frac{c_p^0}{R_u} = - \\frac{a_{k} b_{k}^{2} c_{k} d_{k} \\tau^{2} e^{\\frac{b_{k} \\tau}{T_{c}}}}{T_{c}^{2} \\left(c_{k} + d_{k} e^{\\frac{b_{k} \\tau}{T_{c}}}\\right)^{2}}$$\n", "\n", - "Two important identities:\n", - "$$ \\left(\\frac{1}{\\cosh x}\\right) = \\frac{2e^{-x}}{1+e^{-2x}}$$\n", - "$$ \\left(\\frac{1}{\\sinh x}\\right) = \\frac{2e^{-x}}{1-e^{-2x}}$$\n", - "which yield\n", - "$$ \\left(\\frac{1}{\\cosh x}\\right)^2 = \\frac{4e^{-2x}}{(1+e^{-2x})^2}$$\n", - "$$ \\left(\\frac{1}{\\sinh x}\\right)^2 = \\frac{4e^{-2x}}{(1-e^{-2x})^2}$$" + "Constant is separated out, and handled separately. sinh part can be expanded as\n", + "$$B\\left(\\frac{C/T}{\\sinh(C/T)}\\right)^2 = \\frac{B(-2C/T)^2\\exp(-2C/T)}{(1-\\exp(-2C/T))^2}$$\n", + "where\n", + "$$a_k = B$$\n", + "$$b_k = -2C$$\n", + "$$c_k = 1$$\n", + "$$d_k = -1$$\n", + "\n", + "cosh part can be expanded as\n", + "$$D\\left(\\frac{E/T}{\\cosh(E/T)}\\right)^2 = \\frac{D(-2E/T)^2\\exp(-2E/T)}{(1+\\exp(-2E/T))^2}$$\n", + "where\n", + "$$a_k = -D$$\n", + "$$b_k = -2E$$\n", + "$$c_k = 1$$\n", + "$$d_k = 1$$" ] }, { @@ -102,7 +220,7 @@ "display(diff1.rewrite(exp).subs(B, 2*D).simplify())\n", "term2 = (2*D/T)**2*exp(2*D/T)/(exp(2*D/T)-1)**2\n", "display(term2)\n", - "print 'which is the same as the Aly-Lee sinh term with b_k = 2*D'\n", + "print 'which is the same as the Aly-Lee sinh term with a_k = C, b_k = 2*D, c_k = -1, d_k = 1'\n", "\n", "B, F, T = symbols('B, F, T', real = True, positive = True)\n", "# From Aly-Lee:\n", diff --git a/include/Helmholtz.h b/include/Helmholtz.h index 3d0edd48..7e7cc20f 100644 --- a/include/Helmholtz.h +++ b/include/Helmholtz.h @@ -14,25 +14,24 @@ Residual Helmholtz Energy Terms: Term | Helmholtz Energy Contribution ---------- | ------------------------------ -ResidualHelmholtzPower | \f$ \alpha_r=\left\lbrace\begin{array}{cc}\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} & l_i=0\\ \displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) & l_i\neq 0\end{array}\right.\f$ -ResidualHelmholtzExponential | \f$ \alpha_r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\gamma_i\delta^{l_i}) \f$ -ResidualHelmholtzGaussian | \f$ \alpha_r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\eta_i(\delta-\epsilon_i)^2-\beta_i(\tau-\gamma_i)^2)\f$ -ResidualHelmholtzGERG2008Gaussian | \f$ \alpha_r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\eta_i(\delta-\epsilon_i)^2-\beta_i(\delta-\gamma_i))\f$ -ResidualHelmholtzLemmon2005 | \f$ \alpha_r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) \exp(-\tau^{m_i})\f$ -ResidualHelmholtzNonAnalytic | \f$ \begin{array}{c}\alpha_r&=&\displaystyle\sum_i n_i \Delta^{b_i}\delta\psi \\ \Delta & = & \theta^2+B_i[(\delta-1)^2]^{a_i}\\ \theta & = & (1-\tau)+A_i[(\delta-1)^2]^{1/(2\beta_i)}\\ \psi & = & \exp(-C_i(\delta-1)^2-D_i(\tau-1)^2) \end{array}\f$ -ResidualHelmholtzSAFTAssociating | \f$ \alpha_r = am\left(\ln X-\frac{X}{2}+\frac{1}{2}\right); \f$ +ResidualHelmholtzPower | \f$ \alpha^r=\left\lbrace\begin{array}{cc}\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} & l_i=0\\ \displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) & l_i\neq 0\end{array}\right.\f$ +ResidualHelmholtzExponential | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\gamma_i\delta^{l_i}) \f$ +ResidualHelmholtzGaussian | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\eta_i(\delta-\epsilon_i)^2-\beta_i(\tau-\gamma_i)^2)\f$ +ResidualHelmholtzGERG2008Gaussian | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\eta_i(\delta-\epsilon_i)^2-\beta_i(\delta-\gamma_i))\f$ +ResidualHelmholtzLemmon2005 | \f$ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) \exp(-\tau^{m_i})\f$ +ResidualHelmholtzNonAnalytic | \f$ \begin{array}{c}\alpha^r&=&\displaystyle\sum_i n_i \Delta^{b_i}\delta\psi \\ \Delta & = & \theta^2+B_i[(\delta-1)^2]^{a_i}\\ \theta & = & (1-\tau)+A_i[(\delta-1)^2]^{1/(2\beta_i)}\\ \psi & = & \exp(-C_i(\delta-1)^2-D_i(\tau-1)^2) \end{array}\f$ +ResidualHelmholtzSAFTAssociating | \f$ \alpha^r = am\left(\ln X-\frac{X}{2}+\frac{1}{2}\right); \f$ Ideal-Gas Helmholtz Energy Terms: Term | Helmholtz Energy Contribution ---------- | ------------------------------ -IdealHelmholtzLead | \f$ \alpha_0 = n_1 + n_2\tau + \ln\delta \f$ -IdealHelmholtzEnthalpyEntropyOffset | \f$ \alpha_0 = \displaystyle\frac{\Delta s}{R_u/M}+\frac{\Delta h}{(R_u/M)T}\tau \f$ -IdealHelmholtzLogTau | \f$ \alpha_0 = n_1\log\tau \f$ -IdealHelmholtzPower | \f$ \alpha_0 = \displaystyle\sum_i n_i\tau^{t_i} \f$ -IdealHelmholtzPlanckEinstein | \f$ \alpha_0 = \displaystyle\sum_i n_i\log[1-\exp(-\theta_i\tau)] \f$ -IdealHelmholtzPlanckEinstein2 | \f$ \alpha_0 = \displaystyle\sum_i n_i\log[c_i+\exp(\theta_i\tau)] \f$ +IdealHelmholtzLead | \f$ \alpha^0 = n_1 + n_2\tau + \ln\delta \f$ +IdealHelmholtzEnthalpyEntropyOffset | \f$ \alpha^0 = \displaystyle\frac{\Delta s}{R_u/M}+\frac{\Delta h}{(R_u/M)T}\tau \f$ +IdealHelmholtzLogTau | \f$ \alpha^0 = n_1\log\tau \f$ +IdealHelmholtzPower | \f$ \alpha^0 = \displaystyle\sum_i n_i\tau^{t_i} \f$ +IdealHelmholtzPlanckEinsteinGeneralized | \f$ \alpha^0 = \displaystyle\sum_i n_i\log[c_i+d_i\exp(\theta_i\tau)] \f$ */ class BaseHelmholtzTerm{ public: @@ -107,7 +106,7 @@ struct ResidualHelmholtzPowerElement /*! Term of the form \f[ -\alpha_r=\left\lbrace\begin{array}{cc}\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} & l_i=0\\ \displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) & l_i\neq 0\end{array}\right. +\alpha^r=\left\lbrace\begin{array}{cc}\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} & l_i=0\\ \displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\delta^{l_i}) & l_i\neq 0\end{array}\right. \f] */ class ResidualHelmholtzPower : public BaseHelmholtzTerm{ @@ -162,7 +161,7 @@ struct ResidualHelmholtzExponentialElement }; /** Term of the form -\f[ \alpha_r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\gamma_i\delta^{l_i}) \f] +\f[ \alpha^r=\displaystyle\sum_i n_i \delta^{d_i} \tau^{t_i} \exp(-\gamma_i\delta^{l_i}) \f] */ class ResidualHelmholtzExponential : public BaseHelmholtzTerm{ @@ -573,7 +572,7 @@ public: /// The leading term in the EOS used to set the desired reference state /** \f[ -\alpha_0 = \log(\delta)+a_1+a_2\tau +\alpha^0 = \log(\delta)+a_1+a_2\tau \f] */ class IdealHelmholtzLead : public BaseHelmholtzTerm{ @@ -632,7 +631,7 @@ public: /// The term in the EOS used to shift the reference state of the fluid /** \f[ -\alpha_0 = a_1+a_2\tau +\alpha^0 = a_1+a_2\tau \f] */ class IdealHelmholtzEnthalpyEntropyOffset : public BaseHelmholtzTerm{ @@ -682,7 +681,7 @@ public: /** \f[ -\alpha_0 = a_1\ln\tau +\alpha^0 = a_1\ln\tau \f] */ class IdealHelmholtzLogTau : public BaseHelmholtzTerm @@ -735,7 +734,7 @@ public: /** \f[ -\alpha_0 = \displaystyle\sum_i n_i\tau^{t_i} +\alpha^0 = \displaystyle\sum_i n_i\tau^{t_i} \f] */ class IdealHelmholtzPower : public BaseHelmholtzTerm{ @@ -793,33 +792,42 @@ public: /** \f[ -\alpha_0 = \displaystyle\sum_i n_i\log[1-\exp(-\theta_i\tau)] +\alpha^0 = \displaystyle\sum_i n_i\log[c_i+d_i\exp(\theta_i\tau)] \f] */ -class IdealHelmholtzPlanckEinstein : public BaseHelmholtzTerm{ +class IdealHelmholtzPlanckEinsteinGeneralized : public BaseHelmholtzTerm{ private: - std::vector n,theta; // Use these variables internally + std::vector n,theta,c,d; // Use these variables internally std::size_t N; bool enabled; public: - IdealHelmholtzPlanckEinstein(){N = 0; enabled = false;} + IdealHelmholtzPlanckEinsteinGeneralized(){N = 0; enabled = false;} // Constructor with std::vector instances - IdealHelmholtzPlanckEinstein(std::vector n, std::vector theta) - :n(n), theta(theta) + IdealHelmholtzPlanckEinsteinGeneralized(std::vector n, std::vector theta, std::vector c, std::vector d) + :n(n), theta(theta), c(c), d(d) { N = n.size(); enabled = true; }; - //Destructor - ~IdealHelmholtzPlanckEinstein(){}; + // Destructor + ~IdealHelmholtzPlanckEinsteinGeneralized(){}; + + // Extend the vectors to allow for multiple instances feeding values to this function + void extend(std::vector n, std::vector theta, std::vector c, std::vector d) + { + this->n.insert(this->n.end(), n.begin(), n.end()); + this->theta.insert(this->theta.end(), theta.begin(), theta.end()); + this->c.insert(this->c.end(), c.begin(), c.end()); + this->d.insert(this->d.end(), d.begin(), d.end()); + } bool is_enabled(){return enabled;}; void to_json(rapidjson::Value &el, rapidjson::Document &doc) { - el.AddMember("type","IdealHelmholtzPlanckEinstein",doc.GetAllocator()); + el.AddMember("type","IdealHelmholtzPlanckEinsteinGeneralized",doc.GetAllocator()); cpjson::set_long_double_array("n",n,el,doc); cpjson::set_long_double_array("theta",theta,el,doc); }; @@ -828,81 +836,22 @@ public: long double base(const long double &tau, const long double &delta) throw(){ if (!enabled){return 0.0;} long double s=0; for (std::size_t i=0; i < N; ++i){ - s += n[i]*log(1.0-exp(-theta[i]*tau)); - } return s; + s += n[i]*log(c[i]+d[i]*exp(theta[i]*tau)); + } + return s; }; long double dTau(const long double &tau, const long double &delta) throw(){ if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*theta[i]*(1.0/(1.0-exp(-theta[i]*tau))-1.0);} return s; + long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*theta[i]*d[i]*exp(theta[i]*tau)/(c[i]+d[i]*exp(theta[i]*tau));} + return s; }; long double dTau2(const long double &tau, const long double &delta) throw(){ if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s -= n[i]*pow(theta[i],2)*exp(theta[i]*tau)/pow(1.0-exp(theta[i]*tau),2);} return s; + long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*pow(theta[i],2)*c[i]*d[i]*exp(theta[i]*tau)/pow(c[i]+d[i]*exp(theta[i]*tau),2);} return s; }; long double dTau3(const long double &tau, const long double &delta) throw(){ if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*pow(theta[i],2)*theta[i]*exp(theta[i]*tau)*(exp(theta[i]*tau)+1)/pow(exp(theta[i]*tau)-1,3);} return s; - }; - long double dDelta(const long double &tau, const long double &delta) throw(){return 0.0;}; - long double dDelta2(const long double &tau, const long double &delta) throw(){return 0.0;}; - long double dDelta2_dTau(const long double &tau, const long double &delta) throw(){return 0.0;}; - long double dDelta_dTau(const long double &tau, const long double &delta) throw(){return 0.0;}; - long double dDelta_dTau2(const long double &tau, const long double &delta) throw(){return 0.0;}; - long double dDelta3(const long double &tau, const long double &delta) throw(){return 0;}; -}; - -/** -\f[ -\alpha_0 = \displaystyle\sum_i n_i\log[c_i+\exp(\theta_i\tau)] -\f] -*/ -class IdealHelmholtzPlanckEinstein2 : public BaseHelmholtzTerm{ - -private: - std::vector n,theta,c; // Use these variables internally - std::size_t N; - bool enabled; -public: - IdealHelmholtzPlanckEinstein2(){N = 0; enabled = false;} - // Constructor with std::vector instances - IdealHelmholtzPlanckEinstein2(const std::vector &n, - const std::vector &theta, - const std::vector &c) - :n(n), theta(theta), c(c) - { - N = n.size(); - enabled = true; - }; - - //Destructor - ~IdealHelmholtzPlanckEinstein2(){}; - - bool is_enabled(){return enabled;}; - - void to_json(rapidjson::Value &el, rapidjson::Document &doc) - { - el.AddMember("type","IdealHelmholtzPlanckEinstein2",doc.GetAllocator()); - cpjson::set_long_double_array("n",n,el,doc); - cpjson::set_long_double_array("theta",theta,el,doc); - cpjson::set_long_double_array("c",c,el,doc); - }; - - // Term and its derivatives - long double base(const long double &tau, const long double &delta) throw(){ - if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*log(c[i]+exp(theta[i]*tau));} return s; - }; - long double dTau(const long double &tau, const long double &delta) throw(){ - if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*theta[i]*exp(tau*theta[i])/(c[i]+exp(theta[i]*tau));} return s; - }; - long double dTau2(const long double &tau, const long double &delta) throw(){ - if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*pow(theta[i],2)*c[i]*exp(tau*theta[i])/pow(c[i]+exp(tau*theta[i]),2);} return s; - }; - long double dTau3(const long double &tau, const long double &delta) throw(){ - if (!enabled){return 0.0;} - long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*pow(theta[i],2)*c[i]*(-theta[i])*exp(theta[i]*tau)*(exp(theta[i]*tau)-c[i])/pow(exp(theta[i]*tau)+c[i],3);} return s; + long double s=0; for (std::size_t i=0; i < N; ++i){s += n[i]*pow(theta[i],3)*c[i]*d[i]*(c[i]-d[i]*exp(theta[i]*tau))*exp(theta[i]*tau)/pow(c[i]+d[i]*exp(theta[i]*tau),3);} return s; }; long double dDelta(const long double &tau, const long double &delta) throw(){return 0.0;}; long double dDelta2(const long double &tau, const long double &delta) throw(){return 0.0;}; @@ -986,6 +935,13 @@ public: N = c.size(); }; + void extend(const std::vector &c, const std::vector &t) + { + this->c.insert(this->c.end(), c.begin(), c.end()); + this->t.insert(this->t.end(), t.begin(), t.end()); + N += c.size(); + } + /// Destructor ~IdealHelmholtzCP0PolyT(){}; @@ -1147,100 +1103,89 @@ public: IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffset; IdealHelmholtzLogTau LogTau; IdealHelmholtzPower Power; - IdealHelmholtzPlanckEinstein PlanckEinstein; - IdealHelmholtzPlanckEinstein2 PlanckEinstein2; + IdealHelmholtzPlanckEinsteinGeneralized PlanckEinstein; IdealHelmholtzCP0Constant CP0Constant; IdealHelmholtzCP0PolyT CP0PolyT; - IdealHelmholtzCP0AlyLee CP0AlyLee; long double base(const long double &tau, const long double &delta) { return (Lead.base(tau, delta) + EnthalpyEntropyOffset.base(tau, delta) + LogTau.base(tau, delta) + Power.base(tau, delta) - + PlanckEinstein.base(tau, delta) + PlanckEinstein2.base(tau, delta) + + PlanckEinstein.base(tau, delta) + CP0Constant.base(tau, delta) + CP0PolyT.base(tau, delta) - + CP0AlyLee.base(tau, delta) ); }; long double dDelta(const long double &tau, const long double &delta) { return (Lead.dDelta(tau, delta) + EnthalpyEntropyOffset.dDelta(tau, delta) + LogTau.dDelta(tau, delta) + Power.dDelta(tau, delta) - + PlanckEinstein.dDelta(tau, delta) + PlanckEinstein2.dDelta(tau, delta) + + PlanckEinstein.dDelta(tau, delta) + CP0Constant.dDelta(tau, delta) + CP0PolyT.dDelta(tau, delta) - + CP0AlyLee.dDelta(tau, delta) ); }; long double dTau(const long double &tau, const long double &delta) { return (Lead.dTau(tau, delta) + EnthalpyEntropyOffset.dTau(tau, delta) + LogTau.dTau(tau, delta) + Power.dTau(tau, delta) - + PlanckEinstein.dTau(tau, delta) + PlanckEinstein2.dTau(tau, delta) + + PlanckEinstein.dTau(tau, delta) + CP0Constant.dTau(tau, delta) + CP0PolyT.dTau(tau, delta) - + CP0AlyLee.dTau(tau, delta) ); }; long double dDelta2(const long double &tau, const long double &delta) { return (Lead.dDelta2(tau, delta) + EnthalpyEntropyOffset.dDelta2(tau, delta) + LogTau.dDelta2(tau, delta) + Power.dDelta2(tau, delta) - + PlanckEinstein.dDelta2(tau, delta) + PlanckEinstein2.dDelta2(tau, delta) + + PlanckEinstein.dDelta2(tau, delta) + CP0Constant.dDelta2(tau, delta) + CP0PolyT.dDelta2(tau, delta) - + CP0AlyLee.dDelta2(tau, delta) ); }; long double dDelta_dTau(const long double &tau, const long double &delta) { return (Lead.dDelta_dTau(tau, delta) + EnthalpyEntropyOffset.dDelta_dTau(tau, delta) + LogTau.dDelta_dTau(tau, delta) + Power.dDelta_dTau(tau, delta) - + PlanckEinstein.dDelta_dTau(tau, delta) + PlanckEinstein2.dDelta_dTau(tau, delta) + + PlanckEinstein.dDelta_dTau(tau, delta) + CP0Constant.dDelta_dTau(tau, delta) + CP0PolyT.dDelta_dTau(tau, delta) - + CP0AlyLee.dDelta_dTau(tau, delta) ); }; long double dTau2(const long double &tau, const long double &delta) { return (Lead.dTau2(tau, delta) + EnthalpyEntropyOffset.dTau2(tau, delta) + LogTau.dTau2(tau, delta) + Power.dTau2(tau, delta) - + PlanckEinstein.dTau2(tau, delta) + PlanckEinstein2.dTau2(tau, delta) + + PlanckEinstein.dTau2(tau, delta) + CP0Constant.dTau2(tau, delta) + CP0PolyT.dTau2(tau, delta) - + CP0AlyLee.dTau2(tau, delta)); + ); }; long double dDelta3(const long double &tau, const long double &delta) { return (Lead.dDelta3(tau, delta) + EnthalpyEntropyOffset.dDelta3(tau, delta) + LogTau.dDelta3(tau, delta) + Power.dDelta3(tau, delta) - + PlanckEinstein.dDelta3(tau, delta) + PlanckEinstein2.dDelta3(tau, delta) + + PlanckEinstein.dDelta3(tau, delta) + CP0Constant.dDelta3(tau, delta) + CP0PolyT.dDelta3(tau, delta) - + CP0AlyLee.dDelta3(tau, delta) ); }; long double dDelta2_dTau(const long double &tau, const long double &delta) { return (Lead.dDelta2_dTau(tau, delta) + EnthalpyEntropyOffset.dDelta2_dTau(tau, delta) + LogTau.dDelta2_dTau(tau, delta) + Power.dDelta2_dTau(tau, delta) - + PlanckEinstein.dDelta2_dTau(tau, delta) + PlanckEinstein2.dDelta2_dTau(tau, delta) + + PlanckEinstein.dDelta2_dTau(tau, delta) + CP0Constant.dDelta2_dTau(tau, delta) + CP0PolyT.dDelta2_dTau(tau, delta) - + CP0AlyLee.dDelta2_dTau(tau, delta) ); }; long double dDelta_dTau2(const long double &tau, const long double &delta) { return (Lead.dDelta_dTau2(tau, delta) + EnthalpyEntropyOffset.dDelta_dTau2(tau, delta) + LogTau.dDelta_dTau2(tau, delta) + Power.dDelta_dTau2(tau, delta) - + PlanckEinstein.dDelta_dTau2(tau, delta) + PlanckEinstein2.dDelta_dTau2(tau, delta) + + PlanckEinstein.dDelta_dTau2(tau, delta) + CP0Constant.dDelta_dTau2(tau, delta) + CP0PolyT.dDelta_dTau2(tau, delta) - + CP0AlyLee.dDelta_dTau2(tau, delta) ); }; long double dTau3(const long double &tau, const long double &delta) { return (Lead.dTau3(tau, delta) + EnthalpyEntropyOffset.dTau3(tau, delta) + LogTau.dTau3(tau, delta) + Power.dTau3(tau, delta) - + PlanckEinstein.dTau3(tau, delta) + PlanckEinstein2.dTau3(tau, delta) + + PlanckEinstein.dTau3(tau, delta) + CP0Constant.dTau3(tau, delta) + CP0PolyT.dTau3(tau, delta) - + CP0AlyLee.dTau3(tau, delta) ); }; }; diff --git a/src/Backends/Helmholtz/Fluids/FluidLibrary.h b/src/Backends/Helmholtz/Fluids/FluidLibrary.h index e367369d..f4a3263c 100644 --- a/src/Backends/Helmholtz/Fluids/FluidLibrary.h +++ b/src/Backends/Helmholtz/Fluids/FluidLibrary.h @@ -143,20 +143,36 @@ protected: long double a = cpjson::get_double(contribution,"a"); EOS.alpha0.LogTau = IdealHelmholtzLogTau(a); } + else if (!type.compare("IdealGasHelmholtzPlanckEinsteinGeneralized")) + { + // Retrieve the values + std::vector n = cpjson::get_long_double_array(contribution["n"]); + std::vector t = cpjson::get_long_double_array(contribution["t"]); + + std::vector c = cpjson::get_long_double_array(contribution["c"]); + std::vector d = cpjson::get_long_double_array(contribution["d"]); + if (EOS.alpha0.PlanckEinstein.is_enabled() == true){ + EOS.alpha0.PlanckEinstein.extend(n, t, c, d); + } + else{ + EOS.alpha0.PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized(n, t, c, d); + } + } else if (!type.compare("IdealGasHelmholtzPlanckEinstein")) { - if (EOS.alpha0.PlanckEinstein.is_enabled() == true){throw ValueError("Cannot add ");} + // Retrieve the values std::vector n = cpjson::get_long_double_array(contribution["n"]); std::vector t = cpjson::get_long_double_array(contribution["t"]); - EOS.alpha0.PlanckEinstein = IdealHelmholtzPlanckEinstein(n, t); - } - else if (!type.compare("IdealGasHelmholtzPlanckEinstein2")) - { - if (EOS.alpha0.PlanckEinstein2.is_enabled() == true){throw ValueError("Cannot add");} - std::vector n = cpjson::get_long_double_array(contribution["n"]); - std::vector t = cpjson::get_long_double_array(contribution["t"]); - std::vector c = cpjson::get_long_double_array(contribution["c"]); - EOS.alpha0.PlanckEinstein2 = IdealHelmholtzPlanckEinstein2(n, t, c); + // Flip the sign of theta + for (std::size_t i = 0; i < t.size(); ++i){ t[i] *= -1;} + std::vector c(n.size(), 1); + std::vector d(c.size(), -1); + if (EOS.alpha0.PlanckEinstein.is_enabled() == true){ + EOS.alpha0.PlanckEinstein.extend(n, t, c, d); + } + else{ + EOS.alpha0.PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized(n, t, c, d); + } } else if (!type.compare("IdealGasHelmholtzCP0Constant")) { @@ -177,11 +193,45 @@ protected: } else if (!type.compare("IdealGasHelmholtzCP0AlyLee")) { - if (EOS.alpha0.CP0AlyLee.is_enabled() == true){std::cout << "Cannot add IdealGasHelmholtzCP0AlyLee\n";} - std::vector c = cpjson::get_long_double_array(contribution["c"]); + + std::vector constants = cpjson::get_long_double_array(contribution["c"]); long double Tc = cpjson::get_double(contribution, "Tc"); long double T0 = cpjson::get_double(contribution, "T0"); - EOS.alpha0.CP0AlyLee = IdealHelmholtzCP0AlyLee(c, Tc, T0); + + // Take the constant term if nonzero and set it as a polyT term + if (fabs(constants[0]) > 1e-14){ + std::vector c(1,constants[0]), t(1,0); + if (EOS.alpha0.CP0PolyT.is_enabled() == true){ + EOS.alpha0.CP0PolyT.extend(c,t); + } + else{ + EOS.alpha0.CP0PolyT = IdealHelmholtzCP0PolyT(c, t, Tc, T0); + } + } + + std::vector n, c, d, t; + if (fabs(constants[1]) > 1e-14){ + // sinh term can be converted by setting a_k = C, b_k = 2*D, c_k = -1, d_k = 1 + n.push_back(constants[1]); + t.push_back(-2*constants[2]/Tc); + c.push_back(1); + d.push_back(-1); + } + + if (fabs(constants[3]) > 1e-14){ + // cosh term can be converted by setting a_k = C, b_k = 2*D, c_k = 1, d_k = 1 + n.push_back(-constants[3]); + t.push_back(-2*constants[4]/Tc); + c.push_back(1); + d.push_back(1); + } + + if (EOS.alpha0.PlanckEinstein.is_enabled() == true){ + EOS.alpha0.PlanckEinstein.extend(n, t, c, d); + } + else{ + EOS.alpha0.PlanckEinstein = IdealHelmholtzPlanckEinsteinGeneralized(n, t, c, d); + } } else if (!type.compare("IdealGasHelmholtzEnthalpyEntropyOffset")) { diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index a071c496..a3caf2fe 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -1715,42 +1715,50 @@ long double HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector &mole_fractions, const long double &tau, const long double &delta, const long double &Tr, const long double &rhor) { + long double val; if (is_pure_or_pseudopure) { if (nTau == 0 && nDelta == 0){ - return components[0]->pEOS->base0(tau, delta); + val = components[0]->pEOS->base0(tau, delta); } else if (nTau == 0 && nDelta == 1){ - return components[0]->pEOS->dalpha0_dDelta(tau, delta); + val = components[0]->pEOS->dalpha0_dDelta(tau, delta); } else if (nTau == 1 && nDelta == 0){ - return components[0]->pEOS->dalpha0_dTau(tau, delta); + val = components[0]->pEOS->dalpha0_dTau(tau, delta); } else if (nTau == 0 && nDelta == 2){ - return components[0]->pEOS->d2alpha0_dDelta2(tau, delta); + val = components[0]->pEOS->d2alpha0_dDelta2(tau, delta); } else if (nTau == 1 && nDelta == 1){ - return components[0]->pEOS->d2alpha0_dDelta_dTau(tau, delta); + val = components[0]->pEOS->d2alpha0_dDelta_dTau(tau, delta); } else if (nTau == 2 && nDelta == 0){ - return components[0]->pEOS->d2alpha0_dTau2(tau, delta); + val = components[0]->pEOS->d2alpha0_dTau2(tau, delta); } else if (nTau == 0 && nDelta == 3){ - return components[0]->pEOS->d3alpha0_dDelta3(tau, delta); + val = components[0]->pEOS->d3alpha0_dDelta3(tau, delta); } else if (nTau == 1 && nDelta == 2){ - return components[0]->pEOS->d3alpha0_dDelta2_dTau(tau, delta); + val = components[0]->pEOS->d3alpha0_dDelta2_dTau(tau, delta); } else if (nTau == 2 && nDelta == 1){ - return components[0]->pEOS->d3alpha0_dDelta_dTau2(tau, delta); + val = components[0]->pEOS->d3alpha0_dDelta_dTau2(tau, delta); } else if (nTau == 3 && nDelta == 0){ - return components[0]->pEOS->d3alpha0_dTau3(tau, delta); + val = components[0]->pEOS->d3alpha0_dTau3(tau, delta); } else { throw ValueError(); } + if (!ValidNumber(val)){ + calc_alpha0_deriv_nocache(nTau,nDelta,mole_fractions,tau,delta,Tr,rhor); + throw ValueError(format("calc_alpha0_deriv_nocache returned invalid number with inputs nTau: %d, nDelta: %d", nTau, nDelta)); + } + else{ + return val; + } } else{ // See Table B5, GERG 2008 from Kunz Wagner, JCED, 2012 diff --git a/src/Helmholtz.cpp b/src/Helmholtz.cpp index e2ceb487..23c9bdc8 100644 --- a/src/Helmholtz.cpp +++ b/src/Helmholtz.cpp @@ -1691,10 +1691,11 @@ long double IdealHelmholtzCP0AlyLee::dTau(const long double &tau, const long dou } long double IdealHelmholtzCP0AlyLee::anti_deriv_cp0_tau2(const long double &tau) { - return -c[0]/tau + 2*c[1]*c[2]/(Tc*(exp(-(2*c[2]*tau)/Tc)-1)) + 2*c[3]*c[4]/(Tc*(exp(-(2*c[4]*tau)/Tc)+1)); + return -c[0]/tau + 2*c[1]*c[2]/Tc/(exp(-2*c[2]*tau/Tc)-1) - 2*c[3]*c[4]/Tc*(exp(2*c[4]*tau/Tc)+1); } long double IdealHelmholtzCP0AlyLee::anti_deriv_cp0_tau(const long double &tau) { + long double lnarg =-1 + exp(-2*c[2]*tau/Tc); long double term1 = c[0]*log(tau); long double term2 = 2*c[1]*c[2]*tau/(-Tc + Tc*exp(-2*c[2]*tau/Tc)) + c[1]*log(-1 + exp(-2*c[2]*tau/Tc)) + 2*c[1]*c[2]*tau/Tc; long double term3 = -c[3]*(Tc*exp(2*c[4]*tau/Tc)*log(exp(2*c[4]*tau/Tc) + 1) + Tc*log(exp(2*c[4]*tau/Tc) + 1) - 2*c[4]*tau*exp(2*c[4]*tau/Tc))/(Tc*(exp(2*c[4]*tau/Tc) + 1)); @@ -1716,7 +1717,6 @@ long double IdealHelmholtzCP0AlyLee::dTau3(const long double &tau, const long do /* IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffset; -IdealHelmholtzCP0AlyLee CP0AlyLee; */ #ifdef ENABLE_CATCH @@ -1728,8 +1728,8 @@ class HelmholtzConsistencyFixture public: long double numerical, analytic; - std::tr1::shared_ptr Lead, LogTau, IGPower, PlanckEinstein, PlanckEinstein2, - CP0Constant, CP0PolyT, CP0AlyLee, Gaussian, Lemmon2005, Power, SAFT, NonAnalytic, Exponential, GERG2008; + std::tr1::shared_ptr Lead, LogTau, IGPower, PlanckEinstein, + CP0Constant, CP0PolyT, Gaussian, Lemmon2005, Power, SAFT, NonAnalytic, Exponential, GERG2008; HelmholtzConsistencyFixture(){ Lead.reset(new CoolProp::IdealHelmholtzLead(1,3)); @@ -1739,12 +1739,8 @@ public: IGPower.reset(new CoolProp::IdealHelmholtzPower(n,t)); } { - std::vector n(4,0), t(4,1); n[0] = 0.1; n[2] = 0.5; t[1] = 1; t[2] = 2; t[3] = 2; - PlanckEinstein.reset(new CoolProp::IdealHelmholtzPlanckEinstein(n, t)); - } - { - std::vector n(4,0), t(4,1), c(4,1); n[0] = -0.1; n[2] = 0.1; t[1] = -1; t[2] = -2; t[3] = 2; - PlanckEinstein2.reset(new CoolProp::IdealHelmholtzPlanckEinstein2(n, t, c)); + std::vector n(4,0), t(4,1), c(4,1), d(4,1); n[0] = 0.1; n[2] = 0.5; t[1] = 1; t[2] = 2; t[3] = 2; + PlanckEinstein.reset(new CoolProp::IdealHelmholtzPlanckEinsteinGeneralized(n, t, c, d)); } { long double T0 = 273.15, Tc = 345.857, c = 1.0578, t = 0.33; @@ -1752,11 +1748,6 @@ public: std::vector(1,t), Tc, T0)); } - { - long double T0 = 518.109977174843, Tc = 645.78, c[] = {56.37158920013201, 118.0111016069331, 1792.1, 82.5909330141469, 786.8}; - CP0AlyLee.reset(new CoolProp::IdealHelmholtzCP0AlyLee(std::vector(c, c+sizeof(c)/sizeof(c[0])), - Tc, T0)); - } CP0Constant.reset(new CoolProp::IdealHelmholtzCP0Constant(4/8.314472,300,250)); { long double beta[] = {1.24, 0.821, 15.45, 2.21, 437, 0.743}, @@ -1871,10 +1862,8 @@ public: else if (!t.compare("LogTau")){return LogTau;} else if (!t.compare("IGPower")){return IGPower;} else if (!t.compare("PlanckEinstein")){return PlanckEinstein;} - else if (!t.compare("PlanckEinstein2")){return PlanckEinstein2;} else if (!t.compare("CP0Constant")){return CP0Constant;} else if (!t.compare("CP0PolyT")){return CP0PolyT;} - else if (!t.compare("CP0AlyLee")){return CP0AlyLee;} else if (!t.compare("Gaussian")){return Gaussian;} else if (!t.compare("Lemmon2005")){return Lemmon2005;} diff --git a/src/main.cxx b/src/main.cxx index ebe4d34a..5f232609 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -154,8 +154,8 @@ int main() } if (1) { - double rrr0 = PropsSI("O","T",350,"Q",0,"REFPROP::Water"); - double rrr2 = PropsSI("O","T",350,"Q",0,"Water"); + double rrr0 = PropsSI("C","T",350,"D",1e-13,"REFPROP::MDM"); + double rrr2 = PropsSI("C","T",350,"D",1e-13,"MDM"); double rrr =0 ; } if (1) @@ -171,7 +171,8 @@ int main() } if (1) { - double h1 = PropsSI("H","T",273.15,"Q",0,"CO2"); + double h1 = PropsSI("S","P",101325,"Q",0,"n-Pentane"); + std::string er = get_global_param_string("errstring"); set_reference_stateS("n-Propane","NBP"); double h2 = PropsSI("H","P",101325,"Q",0,"n-Propane"); diff --git a/wrappers/Python/CoolProp5/CoolProp.pyx b/wrappers/Python/CoolProp5/CoolProp.pyx index 401c9c60..be03e060 100644 --- a/wrappers/Python/CoolProp5/CoolProp.pyx +++ b/wrappers/Python/CoolProp5/CoolProp.pyx @@ -72,16 +72,6 @@ cpdef ndarray_or_iterable(object input): # if retval < 0: # raise ValueError('Unable to set reference state') # -# cpdef add_REFPROP_fluid(string_like FluidName): -# """ -# Add a REFPROP fluid to CoolProp internal structure -# -# example:: -# -# add_REFPROP_fluid("REFPROP-PROPANE") -# -# """ -# _add_REFPROP_fluid(FluidName) # # cpdef long get_Fluid_index(string_like Fluid): # """ @@ -116,9 +106,9 @@ cpdef ndarray_or_iterable(object input): cpdef get_global_param_string(string param): return _get_global_param_string(param) -# -# cpdef get_fluid_param_string(string_like fluid, string_like param): -# return _get_fluid_param_string(fluid, param) + +cpdef get_fluid_param_string(string_like fluid, string_like param): + return _get_fluid_param_string(fluid, param) # # cpdef __Props_err1(in1,in2,errstr): # if not len(errstr) == 0: @@ -155,53 +145,53 @@ cpdef PropsSI(in1, in2, in3, in4, in5, in6, in7 = None): else: return _PropsSI(in1, in2, in3, in4, in5, in6, in7) -# cpdef list FluidsList(): -# """ -# Return a list of strings of all fluid names -# -# Returns -# ------- -# FluidsList : list of strings of fluid names -# All the fluids that are included in CoolProp -# -# Notes -# ----- -# -# Here is an example:: -# -# In [0]: from CoolProp.CoolProp import FluidsList -# -# In [1]: FluidsList() -# -# """ -# return _get_global_param_string("FluidsList").split(',') +cpdef list FluidsList(): + """ + Return a list of strings of all fluid names + + Returns + ------- + FluidsList : list of strings of fluid names + All the fluids that are included in CoolProp + + Notes + ----- + + Here is an example:: + + In [0]: from CoolProp.CoolProp import FluidsList + + In [1]: FluidsList() + + """ + return _get_global_param_string("FluidsList").split(',') -# cpdef get_aliases(str Fluid): -# """ -# Return a comma separated string of aliases for the given fluid -# """ -# cdef bytes _Fluid = Fluid.encode('ascii') -# return [F.encode('ascii') for F in (_get_fluid_param_string(_Fluid,'aliases').encode('ascii')).decode('ascii').split(',')] +cpdef get_aliases(str Fluid): + """ + Return a comma separated string of aliases for the given fluid + """ + cdef bytes _Fluid = Fluid.encode('ascii') + return [F.encode('ascii') for F in (_get_fluid_param_string(_Fluid,'aliases').encode('ascii')).decode('ascii').split(',')] # -# cpdef string get_REFPROPname(str Fluid): -# """ -# Return the REFPROP compatible name for the fluid (only useful on windows) -# -# Some fluids do not use the REFPROP name. For instance, -# ammonia is R717, and propane is R290. You can still can still call CoolProp -# using the name ammonia or R717, but REFPROP requires that you use a limited -# subset of names. Therefore, this function that returns the REFPROP compatible -# name. To then use this to call REFPROP, you would do something like:: -# -# In [0]: from CoolProp.CoolProp import get_REFPROPname, Props -# -# In [1]: Fluid = 'REFPROP-' + get_REFPROPname('R290') -# -# In [2]: Props('D', 'T', 300, 'P', 300, Fluid) -# """ -# cdef bytes _Fluid = Fluid.encode('ascii') -# return _get_fluid_param_string(_Fluid,'REFPROP_name') +cpdef string get_REFPROPname(str Fluid): + """ + Return the REFPROP compatible name for the fluid + + Some fluids do not use the REFPROP name. For instance, + ammonia is R717, and propane is R290. You can still can still call CoolProp + using the name ammonia or R717, but REFPROP requires that you use a limited + subset of names. Therefore, this function that returns the REFPROP compatible + name. To then use this to call REFPROP, you would do something like:: + + In [0]: from CoolProp.CoolProp import get_REFPROPname, PropsSI + + In [1]: get_REFPROPname('R290') + + In [2]: PropsSI('D', 'T', 300, 'P', 300, Fluid) + """ + cdef bytes _Fluid = Fluid.encode('ascii') + return _get_fluid_param_string(_Fluid,'REFPROP_name') # cpdef string get_BibTeXKey(str Fluid, str key): # """ @@ -230,31 +220,31 @@ cpdef PropsSI(in1, in2, in3, in4, in5, in6, in7 = None): # Return the current error string # """ # return _get_global_param_string("errstring") -# cpdef set_debug_level(int level): -# """ -# Set the current debug level as integer in the range [0,10] -# -# Parameters -# ---------- -# level : int -# If level is 0, no output will be written to screen, if >0, -# some output will be written to screen. The larger level is, -# the more verbose the output will be -# """ -# _set_debug_level(level) +cpdef set_debug_level(int level): + """ + Set the current debug level as integer in the range [0,10] + + Parameters + ---------- + level : int + If level is 0, no output will be written to screen, if >0, + some output will be written to screen. The larger level is, + the more verbose the output will be + """ + _set_debug_level(level) -# cpdef get_debug_level(): -# """ -# Return the current debug level as integer -# -# Returns -# ------- -# level : int -# If level is 0, no output will be written to screen, if >0, -# some output will be written to screen. The larger level is, -# the more verbose the output will be -# """ -# return _get_debug_level() +cpdef int get_debug_level(): + """ + Return the current debug level as integer + + Returns + ------- + level : int + If level is 0, no output will be written to screen, if >0, + some output will be written to screen. The larger level is, + the more verbose the output will be + """ + return _get_debug_level() # # cpdef bint IsFluidType(string Fluid, string Type): # """ @@ -286,8 +276,8 @@ cdef dict paras = {iDmass : 'D', iCpmass : 'C', # iC0 : 'C0', iCvmass : 'O', -# iV : 'V', -# iL : 'L', + iV : 'V', + iL : 'L', iSmass : 'S', iUmass : 'U', # iDpdT : 'dpdT' @@ -307,16 +297,16 @@ cdef class State: to calculate the enthalpy and pressure. Since the Equations of State are all explicit in temperature and density, each time you call something like:: - h = Props('H','T',T','P',P,Fluid) - s = Props('S','T',T','P',P,Fluid) + h = PropsSI('H','T',T','P',P,Fluid) + s = PropsSI('S','T',T','P',P,Fluid) the solver is used to carry out the T-P flash calculation. And if you wanted entropy as well you could either intermediately calculate ``T``, ``rho`` and then use ``T``, ``rho`` in the EOS in a manner like:: - rho = Props('D','T',T','P',P,Fluid) - h = Props('H','T',T','D',rho,Fluid) - s = Props('S','T',T','D',rho,Fluid) + rho = PropsSI('D','T',T','P',P,Fluid) + h = PropsSI('H','T',T','D',rho,Fluid) + s = PropsSI('S','T',T','D',rho,Fluid) Instead in this class all that is handled internally. So the call to update sets the internal variables in the most computationally efficient way possible @@ -332,7 +322,7 @@ cdef class State: phase : string DEPRECATED : this input is ignored backend : string - The CoolProp backend that should be used + The CoolProp backend that should be used, one of "HEOS" (default), "REFPROP", "INCOMP", "BRINE", etc. """ cdef string _Fluid = Fluid