BioSTEAMDevelopmentGroup / Bioindustrial-Park

BioSTEAM's Premier Repository for Biorefinery Models and Results
MIT License
38 stars 18 forks source link

Error consultation about liquid phase composition #88

Closed zasddsgg closed 11 months ago

zasddsgg commented 1 year ago

Hello, when I called Flash, I encountered the following error, which seems to be caused by denominator being zero. I read the thermosteam.equilibrium.dew_point source code at https://biosteam.readthedocs.io/en/latest/_modules/thermosteam/equilibrium/dew_point.html#:~:text=gamma%20%3D%20f_gamma(x%2C%20T%2C%20*gamma_args)%0A%20%20%20%20denominator%20%3D%20gamma. I found denominator, i.e. gamma, but has not yet found a description of f_gamma. So I don't know how to solve this error. Could I trouble you help to have a look this error? Thanks for your help. Wish you a good day.

Strangely enough, when I execute code "feed.vle(T=180+273.15, P=101000)" instead of Flash, no error is reported.

The code is as follows:

import biosteam as bst
from biosteam import settings
from biosteam import units

Riboflavin = bst.Chemical('Riboflavin', Hf=-55700, Tc=650+273.15, Tb=240+273.15, Hvap=63457)
Ethanol = bst.Chemical('Ethanol')
Water = bst.Chemical('Water')
carbon = bst.Chemical('carbon')
phenol = bst.Chemical('phenol')

bst.settings.set_thermo([Riboflavin, Ethanol, Water, carbon, phenol])

feed = bst.Stream(
    ID='feed',
    price=0.0816,
    units='kg/hr',
    Water=22,
    Ethanol=788,
    carbon=50,
    Riboflavin=30,
    phenol=30,
    T=300+273.15,
    P=7500000,
)

F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), T=180+273.15, P=101000) 
# if use feed.vle(T=180+273.15, P=101000) instead of Flash, no error occurred
F101.simulate()

The error information is as follows:

FloatingPointError                        Traceback (most recent call last)
File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:32, in x_iter(x, x_gamma, T, P, f_gamma, gamma_args)
     31 try:
---> 32     x = x_gamma / denominator
     33 except: 

FloatingPointError: divide by zero encountered in divide

During handling of the above exception, another exception occurred:

Exception                                 Traceback (most recent call last)
Cell In[1], line 27
     13 feed = bst.Stream(
     14     ID='feed',
     15     price=0.0816,
   (...)
     23     P=7500000,
     24 )
     26 F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), T=180+273.15, P=101000)
---> 27 F101.simulate()

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\_unit.py:1732, in Unit.simulate(self, run, design_kwargs, cost_kwargs)
   1730     self._load_stream_links()
   1731     self.run()
-> 1732 self._summary(design_kwargs, cost_kwargs)

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\_unit.py:1465, in Unit._summary(self, design_kwargs, cost_kwargs, lca_kwargs)
   1463 if not (self._design or self._cost): return
   1464 if not self._skip_simulation_when_inlets_are_empty or not all([i.isempty() for i in self._ins]): 
-> 1465     self._design(**design_kwargs) if design_kwargs else self._design()
   1466     self._cost(**cost_kwargs) if cost_kwargs else self._cost()
   1467     self._lca(**lca_kwargs) if lca_kwargs else self._lca()

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\units\_flash.py:291, in Flash._design(self)
    289     self.heat_exchanger._setup() # Removes results
    290 else:
--> 291     self.heat_exchanger.simulate_as_auxiliary_exchanger(self.ins, self.outs, P=self.ins[0].P)

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\units\heat_exchange.py:459, in HXutility.simulate_as_auxiliary_exchanger(self, ins, outs, duty, vle, scale, hxn_ok, P)
    457 inlet.mix_from(ins, energy_balance=False)
    458 if P is None: P = inlet.P
--> 459 if vle: inlet.vle(H=sum([i.H for i in ins]), P=P)
    460 if outs is None:
    461     if duty is None: raise ValueError('must pass duty when no outlets are given')

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\vle.py:377, in VLE.__call__(self, T, P, V, H, S, x, y)
    375 elif H_spec:
    376     try:
--> 377         self.set_PH(P, H, stacklevel=1)
    378     except NoEquilibrium:
    379         thermal_condition = self._thermal_condition

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\vle.py:1067, in VLE.set_PH(self, P, H, stacklevel)
   1064     return
   1066 # Check if super heated vapor
-> 1067 T_dew, x_dew = self._dew_point.solve_Tx(self._z, P)
   1068 if T_dew <= T_bubble: 
   1069     T_dew, T_bubble = T_bubble, T_dew

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:228, in DewPoint.solve_Tx(self, z, P)
    226 args = (P, z_norm, zP, x)
    227 try:
--> 228     T = flx.aitken_secant(f, T_guess, T_guess + 1e-3,
    229                           self.T_tol, 5e-12, args,
    230                           checkiter=False)
    231 except RuntimeError:
    232     Tmin = self.Tmin

File D:\anaconda\envs\zddd\lib\site-packages\flexsolve\open_solvers.py:68, in aitken_secant(f, x0, x1, xtol, ytol, args, maxiter, checkroot, checkiter)
     66     dx = -2 * xtol
     67 abs_ = abs
---> 68 y0 = f(x0, *args)
     69 if abs_(y0) < ytol: return x0
     70 aitken_iter = utils.scalar_aitken_iter

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:138, in DewPoint._T_error(self, T, P, z_norm, zP, x)
    136 pcf = self.pcf(T, P, Psats)
    137 x_gamma = phi * zP / Psats / pcf
--> 138 x[:] = self._solve_x(x_gamma, T, P, x)
    139 return 1 - x.sum()

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:129, in DewPoint._solve_x(self, x_gamma, T, P, x)
    127 def _solve_x(self, x_gamma, T, P, x):
    128     gamma = self.gamma
--> 129     return solve_x(x, x_gamma, T, P, gamma.f, gamma.args)

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:45, in solve_x(x_guess, x_gamma, T, P, f_gamma, gamma_args)
     43 def solve_x(x_guess, x_gamma, T, P, f_gamma, gamma_args):
     44     args = (x_gamma, T, P, f_gamma, gamma_args)
---> 45     x = flx.wegstein(x_iter, x_guess, 1e-10, args=args, checkiter=False,
     46                      checkconvergence=False, convergenceiter=3)
     47     return x

File D:\anaconda\envs\zddd\lib\site-packages\flexsolve\iterative_solvers.py:58, in wegstein(f, x, xtol, args, maxiter, checkiter, checkconvergence, convergenceiter)
     56 errors = np.zeros(convergenceiter)
     57 x0 = x
---> 58 x1 = g0 = f(x0, *args)
     59 wegstein_iter = utils.wegstein_iter
     60 fixedpoint_converged = utils.fixedpoint_converged

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:34, in x_iter(x, x_gamma, T, P, f_gamma, gamma_args)
     32     x = x_gamma / denominator
     33 except: 
---> 34     raise Exception('liquid phase composition is infeasible')
     35 if (x < 0).any():
     36     raise Exception('liquid phase composition is infeasible')

Exception: liquid phase composition is infeasible
zasddsgg commented 1 year ago

It seems to be the error caused by the activity coefficient of 0. I tried to change the operating conditions of Flash, such as setting "V=0.02, P=101000", "Q=0, P=101000", but the above error was still reported.

zasddsgg commented 1 year ago

Hello, I tried different conditions for solving the above error, but error was still reported under the conditions I set. If I remove phenol from the feed stream, then no error was reported. But actually, phenol is in my system. It is a little strange. If I do not remove phenol, execute feed.vle(T=180+273.15, P=101000) instead of F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), T=180+273.15, P=101000), there was no error, which was strange. Isn't Flash() the same as Stream.vle()? Why does one report an error and the other not? Does Stream.vle() and Flash(), set the same parameters, such as V,P or T,P or V,T, return the same result? Thanks for your help. Wish you a good day.

zasddsgg commented 12 months ago

Hello, I still haven't solved the above error, could you please have a look at it when you are free. Thanks for your help. Wish you a good day.

yalinli2 commented 11 months ago

from what I can tell the problem is not in the stream's vle, it's in the heat exchanger of the Flash, the F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), Q=0, P=101000) could work, but I'm not sure if it makes any sense for your system.

also, removing phenol from the system won't solve the problem for me (i.e., the following code still trigger the error):

import biosteam as bst
from biosteam import settings
from biosteam import units

Riboflavin = bst.Chemical('Riboflavin', Hf=-55700, Tc=650+273.15, Tb=240+273.15, Hvap=63457)
Ethanol = bst.Chemical('Ethanol')
Water = bst.Chemical('Water')
carbon = bst.Chemical('carbon')
phenol = bst.Chemical('phenol')

bst.settings.set_thermo([Riboflavin, Ethanol, Water, carbon, phenol])

feed = bst.Stream(
    ID='feed',
    price=0.0816,
    units='kg/hr',
    Water=22,
    Ethanol=788,
    carbon=50,
    Riboflavin=30,
    # phenol=30,
    T=300+273.15,
    P=7500000,
)
F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), T=180+273.15, P=7500000) 
zasddsgg commented 11 months ago

Thanks for your reply, I tried using F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), Q=0, P=101000), but the following error was still triggered. Moreover, the conditions for Flash in the process is actually T=180+273.15, P=101000. May I ask you how can solve errors in the case T=180+273.15, P=101000 and Q=0, P=101000? Thanks for your help. Wish you a good day.

The error using F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), Q=0, P=101000) is as follows:

FloatingPointError                        Traceback (most recent call last)
File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:32, in x_iter(x, x_gamma, T, P, f_gamma, gamma_args)
     31 try:
---> 32     x = x_gamma / denominator
     33 except: 

FloatingPointError: divide by zero encountered in divide

During handling of the above exception, another exception occurred:

Exception                                 Traceback (most recent call last)
Cell In[1], line 28
     26 F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), Q=0, P=101000) 
     27 # if use feed.vle(T=180+273.15, P=101000) instead of Flash, no error occurred
---> 28 F101.simulate()

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\_unit.py:1731, in Unit.simulate(self, run, design_kwargs, cost_kwargs)
   1729     for ps in self._specifications: ps.compile_path(self)
   1730     self._load_stream_links()
-> 1731     self.run()
   1732 self._summary(design_kwargs, cost_kwargs)

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\_unit.py:1392, in Unit.run(self)
   1390         if self.run_after_specifications: self._run()
   1391 else:
-> 1392     self._run()

File D:\anaconda\envs\zddd\lib\site-packages\biosteam\units\_flash.py:267, in Flash._run(self)
    266 def _run(self):
--> 267     separations.vle(self.ins[0], *self.outs, self.T, self.P, self.V, 
    268                     self.Q, self.x, self.y, self._multi_stream)

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\separations.py:727, in vle(feed, vap, liq, T, P, V, Q, x, y, multi_stream)
    725     ms = feed.copy()
    726 H = feed.H + Q if Q is not None else None
--> 727 ms.vle(P=P, H=H, T=T, V=V, x=x, y=y)
    729 # Set Values
    730 vap.phase = 'g'

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\vle.py:377, in VLE.__call__(self, T, P, V, H, S, x, y)
    375 elif H_spec:
    376     try:
--> 377         self.set_PH(P, H, stacklevel=1)
    378     except NoEquilibrium:
    379         thermal_condition = self._thermal_condition

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\vle.py:1067, in VLE.set_PH(self, P, H, stacklevel)
   1064     return
   1066 # Check if super heated vapor
-> 1067 T_dew, x_dew = self._dew_point.solve_Tx(self._z, P)
   1068 if T_dew <= T_bubble: 
   1069     T_dew, T_bubble = T_bubble, T_dew

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:228, in DewPoint.solve_Tx(self, z, P)
    226 args = (P, z_norm, zP, x)
    227 try:
--> 228     T = flx.aitken_secant(f, T_guess, T_guess + 1e-3,
    229                           self.T_tol, 5e-12, args,
    230                           checkiter=False)
    231 except RuntimeError:
    232     Tmin = self.Tmin

File D:\anaconda\envs\zddd\lib\site-packages\flexsolve\open_solvers.py:68, in aitken_secant(f, x0, x1, xtol, ytol, args, maxiter, checkroot, checkiter)
     66     dx = -2 * xtol
     67 abs_ = abs
---> 68 y0 = f(x0, *args)
     69 if abs_(y0) < ytol: return x0
     70 aitken_iter = utils.scalar_aitken_iter

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:138, in DewPoint._T_error(self, T, P, z_norm, zP, x)
    136 pcf = self.pcf(T, P, Psats)
    137 x_gamma = phi * zP / Psats / pcf
--> 138 x[:] = self._solve_x(x_gamma, T, P, x)
    139 return 1 - x.sum()

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:129, in DewPoint._solve_x(self, x_gamma, T, P, x)
    127 def _solve_x(self, x_gamma, T, P, x):
    128     gamma = self.gamma
--> 129     return solve_x(x, x_gamma, T, P, gamma.f, gamma.args)

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:45, in solve_x(x_guess, x_gamma, T, P, f_gamma, gamma_args)
     43 def solve_x(x_guess, x_gamma, T, P, f_gamma, gamma_args):
     44     args = (x_gamma, T, P, f_gamma, gamma_args)
---> 45     x = flx.wegstein(x_iter, x_guess, 1e-10, args=args, checkiter=False,
     46                      checkconvergence=False, convergenceiter=3)
     47     return x

File D:\anaconda\envs\zddd\lib\site-packages\flexsolve\iterative_solvers.py:58, in wegstein(f, x, xtol, args, maxiter, checkiter, checkconvergence, convergenceiter)
     56 errors = np.zeros(convergenceiter)
     57 x0 = x
---> 58 x1 = g0 = f(x0, *args)
     59 wegstein_iter = utils.wegstein_iter
     60 fixedpoint_converged = utils.fixedpoint_converged

File D:\anaconda\envs\zddd\lib\site-packages\thermosteam\equilibrium\dew_point.py:34, in x_iter(x, x_gamma, T, P, f_gamma, gamma_args)
     32     x = x_gamma / denominator
     33 except: 
---> 34     raise Exception('liquid phase composition is infeasible')
     35 if (x < 0).any():
     36     raise Exception('liquid phase composition is infeasible')

Exception: liquid phase composition is infeasible
yalinli2 commented 11 months ago

I'm not sure of your environment, but using the master branch of biosteam/thermosteam, if I run

import biosteam as bst
from biosteam import settings
from biosteam import units

Riboflavin = bst.Chemical('Riboflavin', Hf=-55700, Tc=650+273.15, Tb=240+273.15, Hvap=63457)
Ethanol = bst.Chemical('Ethanol')
Water = bst.Chemical('Water')
carbon = bst.Chemical('carbon')
phenol = bst.Chemical('phenol')

bst.settings.set_thermo([Riboflavin, Ethanol, Water, carbon, phenol])

feed = bst.Stream(
    ID='feed',
    price=0.0816,
    units='kg/hr',
    Water=22,
    Ethanol=788,
    carbon=50,
    Riboflavin=30,
    phenol=30,
    T=300+273.15,
    P=7500000,
)

F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), Q=0, P=101000) 
# if use feed.vle(T=180+273.15, P=101000) instead of Flash, no error occurred
F101.simulate()
F101.show()

I get

Flash: F101
ins...
[0] feed  
    phase: 'l', T: 573.15 K, P: 7.5e+06 Pa
    flow (kmol/hr): Riboflavin  0.0797
                    Ethanol     17.1
                    Water       1.22
                    carbon      4.16
                    phenol      0.319
outs...
[0] vapor  
    phase: 'g', T: 1431.7 K, P: 101000 Pa
    flow (kmol/hr): Riboflavin  0.0797
                    Ethanol     17.1
                    Water       1.22
                    carbon      4.16
                    phenol      0.319
[1] liquid  
    phase: 'l', T: 1431.7 K, P: 101000 Pa
    flow: 0

The script run, but design-wise I'm not sure if it makes a lot of sense. You have that carbon chemical with no designation of properties, which could cause problems in simulation.

Also, biosteam's default heating agents are all steams and their temperature limits are not that high, won't be sufficient for your feed stream that's at 573.15 K:

>>> bst.HeatUtility.heating_agents

[<UtilityAgent: low_pressure_steam>,
 <UtilityAgent: medium_pressure_steam>,
 <UtilityAgent: high_pressure_steam>]

>>> for i in bst.HeatUtility.heating_agents: print(i.T)
412.189
454.77
508.991

If you take out that carbon chemical and adjust down your feed temperature, you'll get more reasonable results:

import biosteam as bst
from biosteam import settings
from biosteam import units

Riboflavin = bst.Chemical('Riboflavin', Hf=-55700, Tc=650+273.15, Tb=240+273.15, Hvap=63457)
Ethanol = bst.Chemical('Ethanol')
Water = bst.Chemical('Water')
carbon = bst.Chemical('carbon')
phenol = bst.Chemical('phenol')

bst.settings.set_thermo([Riboflavin, Ethanol, Water, carbon, phenol])

feed = bst.Stream(
    ID='feed',
    price=0.0816,
    units='kg/hr',
    Water=22,
    Ethanol=788,
    # carbon=50,
    Riboflavin=30,
    phenol=30,
    T=220+273.15,
    P=7500000,
)

F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), T=180+273.15, P=101000)
# if use feed.vle(T=180+273.15, P=101000) instead of Flash, no error occurred
F101.simulate()
F101.show()
Flash: F101
ins...
[0] feed  
    phase: 'l', T: 493.15 K, P: 7.5e+06 Pa
    flow (kmol/hr): Riboflavin  0.0797
                    Ethanol     17.1
                    Water       1.22
                    phenol      0.319
outs...
[0] vapor  
    phase: 'g', T: 453.15 K, P: 101000 Pa
    flow (kmol/hr): Riboflavin  0.0797
                    Ethanol     17.1
                    Water       1.22
                    phenol      0.319
[1] liquid  
    phase: 'l', T: 453.15 K, P: 101000 Pa
    flow: 0

But if your feed is indeed at that high temperature, you'll want to add heating agent that can transfer heat at higher temperature

yoelcortes commented 11 months ago

@yalinli2, @zasddsgg thanks for doing detective work to characterize the error! I was able to figure out the following workaround thanks to it!

I'm not sure what models are in carbon but I don't think it makes sense for it to be in liquid or vapor states. Based on the discussion, we know it is messing up your VLE. If you force it to be a solid only (and therefore not participate in phase equilibrium), then the problem is resolved.

import biosteam as bst
from biosteam import settings
from biosteam import units

Riboflavin = bst.Chemical('Riboflavin', Hf=-55700, Tc=650+273.15, Tb=240+273.15, Hvap=63457)
Ethanol = bst.Chemical('Ethanol')
Water = bst.Chemical('Water')
carbon = bst.Chemical('carbon', phase='s') # !!! Make sure carbon does not participate in VLE
phenol = bst.Chemical('phenol')

bst.settings.set_thermo([Riboflavin, Ethanol, Water, carbon, phenol])

feed = bst.Stream(
    ID='feed',
    price=0.0816,
    units='kg/hr',
    Water=22,
    Ethanol=788,
    carbon=50,
    Riboflavin=30,
    phenol=30,
    T=300+273.15,
    P=7500000,
)

F101 = units.Flash('F101', feed, outs=('vapor', 'liquid'), T=180+273.15, P=101000) 
F101.simulate()