the flashTQ function went wrong in calculating Binary mixture: methane-benzene's VLE

[funnypoem]

in test_cython.py line 616 when x0=0.0252, I can get a precise result, but I tried other X(methane) from the source literature such as T=421.05K(147.9℃) , X(methane)=0.0502, 0.0823 and other T, X(methane) as 188.7℃ 228.0℃. Anyway, only minimum X(methane) will return a correct result. In conclusion, for higher Patm in any T, the code would go wrong

source: H.-M. Lin, H. M. Sebastian, J. J. Simnick, and K.-C. Chao, “Gas-liquid equilibrium in binary mixtures of methane with N-decane, benzene, and toluene,” J. Chem. Eng. Data, vol. 24, no. 2, pp. 146–149, Apr. 1979.

[zmeri]

Thank you for reporting this bug. I'll take a look and see if I can figure out what's causing this behavior.

On Fri, Apr 15, 2022 at 6:34 PM funnypoem @.***> wrote:

in test_cython.py line 616 when x0 http://methane=0.0252, I can get a precise result, but I tried other X(methane) from the source literature such as T=421.05K(147.9℃) , X(methane)=0.0502, 0.0823 and other T, X(methane) as 188.7℃ 228.0℃. Anyway, only minimum X(methane) will return a correct result. In conclusion, for higher Patm in any T, the code would go wrong

source: H.-M. Lin, H. M. Sebastian, J. J. Simnick, and K.-C. Chao, “Gas-liquid equilibrium in binary mixtures of methane with N-decane, benzene, and toluene,” J. Chem. Eng. Data, vol. 24, no. 2, pp. 146–149, Apr. 1979.

— Reply to this email directly, view it on GitHub https://github.com/zmeri/PC-SAFT/issues/98, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFKXTOG6MOOUTFGKIFG454TVFGEBZANCNFSM5TQWLFBQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[funnypoem]

I reviewed the cpp code and found that bug might appear in scan_pressure function(pcsaft_electrolyte.cpp). In the for cycle, when "i" increased to 3, the cycle would break, which means "err_min" always less than 1e9 and the "if break" sentence run. If I delete the limit of “ctr_increasing > 2”, the flashTQ function would return a result without any expections. But the result is not correct, I think the reason is that there are more than one volume roots and the code chose a wrong one without the limit of “ctr_increasing > 2”. I am trying to solve this problem. By the way, I'm a I am a PhD student and I've been working on some applications of the SAFTs model recently. I was wondering if you can give me your email address so that I can give you some ideas of my calculations. My email address is songyc97@outlook.com. Thank you!

[funnypoem]

For this flashTQ calculation bug, I have a guess about the code which is located in outerTQ. When estimate_flash_p gave a pguess(usually similar with the target p) that had a higher value(more than 1e6), the outerTQ function would calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.

 rhol = pcsaft_den_cpp(t, p, x, 0, cppargs);
 rhov = pcsaft_den_cpp(t, p, x, 1, cppargs);
 if ((rhol - rhov) < 1e-4) {
     throw SolutionError("liquid and vapor densities are the same.");
 }

Here rhol and rhov used the same x. Although these sentences used different phase number(0 or 1), pcsaft_den_cpp would return the same result and then, rhol would equal to rhov and the code: throw SolutionError("liquid and vapor densities are the same.") run. So, should the same x be avoided here? For example, we can give xl and assume that xv[i]=1/ncomp. I tried this and the calculation got a valid result. I am going to try different literature values next to verify the effectiveness of this approach.

[zmeri]

Yes, I think you are on the right track. I think the problem here is that, as you noticed, methane is so much more volatile than benzene. This means that only a single phase can be found by the density solver when using the overall system composition (x). Although you were able to get an answer by simply setting x to be 1/ncomp, this probably isn't a robust solution, so yes, we need to look for a better approach. Another option, when the components have drastically different volatilities, might be to estimate an initial value for the liquid-vapor distribution ratio (k) using the vapor pressure for each compound at the given temperature (k=Psat/P). However, for the current system the problem is that methane is supercritical and does not have a pure component vapor pressure at the given conditions. I think for cases like this, when gases are dissolved in a liquid at small concentrations, Henry's law is generally used. Maybe that would be a good starting point for an initial guess if the first approach (simply using x) doesn't work. I'm trying a few of these different approaches as well, but I haven't found one yet that I think is good enough. I'll let you know when I make some progress.

On Tue, Apr 19, 2022 at 8:31 PM funnypoem @.***> wrote:

For this flashTQ calculation bug, I have a guess about the code which is located in outerTQ. When estimate_flash_p gave a pguess(usually similar with the target p) that had a higher value(more than 1e6), the outerTQ function would calculate initial guess for compositions based on fugacity coefficients and Raoult's Law.

rhol = pcsaft_den_cpp(t, p, x, 0, cppargs); rhov = pcsaft_den_cpp(t, p, x, 1, cppargs); if ((rhol - rhov) < 1e-4) { throw SolutionError("liquid and vapor densities are the same."); }

Here rhol and rhov used the same x. Although these sentences used different phase number(0 or 1), pcsaft_den_cpp would return the same result and then, rhol would equal to rhov and the code: throw SolutionError("liquid and vapor densities are the same.") run. So, should the same x be avoided here? For example, we can give xl and assume that xv[i]=1/ncomp. I tried this and the calculation got a valid result. I am going to try different literature values next to verify the effectiveness of this approach.

— Reply to this email directly, view it on GitHub https://github.com/zmeri/PC-SAFT/issues/98#issuecomment-1102908813, or unsubscribe https://github.com/notifications/unsubscribe-auth/AFKXTOECOEOYJGOTOS6DYNLVF3UVFANCNFSM5TQWLFBQ . You are receiving this because you commented.Message ID: @.***>

[zmeri]

Ok, it seems I was able to solve the problem. I updated the flash function by adding an if statement that attempts to calculate an initial k value a different way if the default approach is unsuccessful. I still need to clean up the code and finish testing and then I will push an update.

[funnypoem]

Congratulations! I'm glad to be able to solve this problem with you. I'm using PC-SAFT model to do some refrigerant property calculations and trying to develop SAFT-VR/Mie projects based on your code frame. Looking forward more cooperation with you on SAFTs projects in the future.

[zmeri]

I have updated the code with the improvements to the flash algorithm. I added an additional test as well for methane + benzene at a higher mole fraction of methane. Let me know if this resolves the issue you've been having.

However, only the flashTQ algorithm is able to reliably solve the methane + benzene system. Interestingly, the flashPQ algorithm cannot. It appears that the PC-SAFT equation of state exhibits some anomalous behavior for the methane + benzene system. This can be seen when looking at the k ratio for methane in the mixture (see the figure below). For a fixed pressure, as the temperature rises, the k value initially increases. This k ratio represents how readily a compound moves to the vapor phase and it is logical that this ratio would increase with temperature. However, for the methane + benzene system PC-SAFT shows that as successively higher temperatures the k value peaks and then starts to decline.

This doesn't seem realistic to me because why should a compound be less likely to vaporize at higher temperatures. The k values for benzene do not appear to behave this way. Indeed, for most systems I've looked at the k values continue to increase with temperature. The only systems I've seen that show this abnormal maximum and subsequent decrease in k are those with light gases, for example, methane + benzene and hydrogen + hexadecane. It may be that this behavior comes directly from the PC-SAFT equations themselves. It is also possible that the behavior comes from some sort of bug in this specific code, but I've also checked my values for hydrogen + hexadecane with those given by the original PC-SAFT code given by Gross and the results matched, which makes me wonder if this behavior is somehow inherent to the PC-SAFT equations themselves. Privat et al. have shown before that the PC-SAFT equations can result in more than 3 density roots and can lead to spurious phase equilibrium points. Maybe this behavior is related.

The flash algorithm I am using was proposed by Watson et al. (https://doi.org/10.1021/acs.iecr.6b03956), and this algorithm essentially assumes that the equilibrium ratio between the phases increases linearly with 1/T. The algorithm uses this kind of a linear equation to select a new temperature during each iteration. From what I have seen, it seems that this flash algorithm cannot find the correct flash temperature for systems like methane + benzene because the k value for methane starts decreasing at some point. Thus, the flashPQ algorithm fails. Interestingly, the k value for methane still appears to continually increase with pressure, and this may be why the flashTQ algorithm still succeeds.

To get the flashPQ algorithm to also work for the methane + benzene system, a lot more work would need to be done to understand why this anomalous behavior occurs and see if there is a way to adjust the flash algorithm to allow it to still find the correct solution.

[Nram94]

I have been following these exchanges and regarding your last comment I have a question: how are you calculating the k_values for methane? The y-axis title indicates y/x =fug_coef_liq/fug_coef_vap is this correct?

on the other hand it may be that you are on a fictious saturation line for methane as noted by Privat et al. If this is the issue the soultion is not as straigthforward.

However I have noticed that you are working with temperatures well above the critical temperature of methane. The PC-SAFT struggles inherently to calculate the critical pressures (deviations can be around 20% on average). It may be that the supercritical properties are also inaccurately calculated.

[Nram94]

I just checked regarding just methane the PC-SAFT critical pressure deviation is around 17%.

Here is the reference: : 10.1021/acs.jced.0c00792

[zmeri]

@Nram94 I am calculating the k values this way: k = fugcoef_l / fugcoef_v, as shown in the y-axis label in the figure.

I don't think that the issue is caused by a fictitious saturation line for methane because from what I remember of Privat et al.'s article the fictitious saturation line occurs usually at lower temperatures. Although, I may be wrong because I haven't checked this to be sure.

And yes, as you mentioned, I suspect as well that the problem might be related to the fact that in the mixture methane is well above its critical temperature. I have also heard before that PC-SAFT doesn't perform as well in the supercritical region.

[Nram94]

@zmeri I am inclined to think that the issue is the latter. One way to test this is to parametrize methane (and perhaps benzene too) by exact reproduction of the critical point and the acentric factor. The publicarion I referenced has the parametwrs for both methane and benzene doing this parametrization (with the addition of a volume translation). This publication has more details on the procedure : https://doi.org/10.1021/acs.iecr.9b04660

I should mention that if you only reproduce Tc, Pc and omega without volume translation, the accuracy on the liquid density suffers notably (sort of like the cubic EoS do). I do not have time at this moment to do the test myself, but it may give some insight on whether the problem is the model or something else.

[funnypoem]

The nonlinearity of equilibrium ratio k with 1/T might cause flashPQ algorithm cannot solve the system reliably, The limitations of the model in the supercritical region may need to be evaluated. However, the flash algorithm has other problems and I had confirmed one in calculation. I noticed that flashTQ calculations for a particular binary system would return the same xl and xv with a unreliable P. I tested the CO2 - n-decane system at 238C from the Gross's publication(10.1021/ie0003887) We can see that the algorithm returned the same xl and xv in high pressure and it's hard to avoid the bug. In the Iterative calculation, when the algorithm give a similar xl and xv in any pressure, ki would close to 1 and ln(ki) close to 0. As a result, lnkb is close to 0, In this case, both A and B are equal to 0, and the value of pressure is irrelevant because the iterative will converge. I initially guessed that this convergence error was related to the initial pressure value p_guess of outerTQ, as estimate_flash_p algorithm and pressure traverse from 10^-6 ~ 10^9 give a unreliable p_guess. But when I tested the whole pressure range, I found that all of the p_guess in pressure traverse from 10^-6 ~ 10^9 will return the same result. I try to calculate with Gross's fortran program and also have the same problem.

[zmeri]

I am seeing similar behavior for some points, but for others I get reasonable results. Here is my curve.

When I checked one of the points where xl and xv were the same, I was able to get the correct result by supplying an accurate guess value to the flashTQ function (the results in the figure above are without supplying an initial guess). Now the flashTQ function should ideally be able to find the correct result without needing the user to supply a guess, so I agree that the algorithm could be improved. This does at least indicate, though, that the problem is in how the initial guess is calculated. It might be a problem with finding an initial p_guess, but it seems more likely to me that the issue is with finding initial ki values for each of the components.

I will try to take a look at this and see if I can find a better way to initialize the k values, but I must admit that I have a feeling this will take a fair amount of work. I don't have a lot of time to devote to this since this is just a side project, so it may be a while before I can find a better solution. If you make some progress, let me know. Maybe you'll find a solution before I do.

[gladio145]

Is this issue fixed? @zmeri

[zmeri]

This doesn't require a fix, it requires a new feature or potentially changes to PC-SAFT itself. The core PC-SAFT functions work correctly, as in they match with the results of other implementations and literature data as far as I've tested. One problem, as mentioned here, was that the k value for light gases in these mixtures (methane, hydrogen) does not increase linearly with temperature, which causes the PQ flash to fail. The TQ flash can still work, but sometimes fails if the initial guess is not close enough. I played around a little with different ways to get an initial guess, but didn't find a robust way to get a better guess.

These issues seem to require thorough studying to get good results for these sorts of mixtures with light gases. Or, potentially, a different flash algorithm would be able to perform better. However, these solutions would require significant development and I do not have the time to do that work. If someone else is willing to expand the code, I can incorporate a pull request with the improvements. Another option would be to try other PC-SAFT implementations that use different flash algorithms, such as Clapeyron or teqp. I don't think those implementations currently include association, polar, or electrolyte terms, but for mixtures like methane and benzene those other terms aren't necessary.

[PPG2889]

in test_cython.py line 616 when x0=0.0252, I can get a precise result, but I tried other X(methane) from the source literature such as T=421.05K(147.9℃) , X(methane)=0.0502, 0.0823 and other T, X(methane) as 188.7℃ 228.0℃. Anyway, only minimum X(methane) will return a correct result. In conclusion, for higher Patm in any T, the code would go wrong

source: H.-M. Lin, H. M. Sebastian, J. J. Simnick, and K.-C. Chao, “Gas-liquid equilibrium in binary mixtures of methane with N-decane, benzene, and toluene,” J. Chem. Eng. Data, vol. 24, no. 2, pp. 146–149, Apr. 1979.

I'm also facing similar issue. I've attached code with output here.

[PPG2889]

Ok, it seems I was able to solve the problem. I updated the flash function by adding an if statement that attempts to calculate an initial k value a different way if the default approach is unsuccessful. I still need to clean up the code and finish testing and then I will push an update.

Hi, I'm also facing the similar issue as reported by @funnypoem. I've attached here screenshot of my program and output. please have a look.

Regards, Pradnya

[zmeri]

@PPG2889 This problem has already been discussed in previous comments on this issue. As mentioned, mixtures with light gases are difficult to perform flash calculations for and most likely the problem is with the way the flash algorithm calculates an initial guess. You could test this by providing the flashTQ function with an initial guess for the pressure. Solving this problem would require significant development and I am not able to contribute that time right now. If you would like to contribute an improved flash algorithm, I can incorporate it into the code here.