Equilibrium calculation for ammonia - H2O flash

[bleves91]

Hello I have the following code to compute bubble temp. at pressure P, with assigned composition:

from iapws import H2ONH3
from iapws import ammonia as iapws
from scipy.optimize import fsolve

mix=iapws.H2ONH3()

#----------------- #fugacity in funciton of ptx
def fugNH3(p, t, x):
def hptx(rho):
y=mix._prop(rho, t, x)["P"]-p
return y
rho=fsolve(hptx, 500)[0]
return mix._prop(rho, t, x)["fugNH3"]

#----------------- #bubble temp
def bubble(p, z, tinit, y1init):

x1=z
x2=z-1
y1=y1init
t=tinit
res=1
while res>0.00001:
    y2=1-y1
    k1=fugNH3(p, t, x1)/fugNH3(p, t, y1)
    k2=(1-k1*x1)/x2
    y1=k1*x1
    tinit=t
    def eq(t):  
        eqn=fugNH3(p, t, x1)/fugNH3(p, t, y1)*x1+(1-k1*x1)/x2*x2-1
        return eqn
    t=fsolve(eq, tinit)[0]
    res=t-tinit
return t 

print(bubble(3.5, 0.5, 110+273, 0.95))

I am struggling to get any results. I wonder, the property "fugNH3" is the fugacity or the fugacity coefficient?
Please help me if you have faced some kind of problem like this

[jjgomera]

Hi, excuse me for my late response

Try with this code

from math import log, exp
from numpy import arange
from scipy.optimize import fsolve
import iapws
from iapws.ammonia import NH3


mix = iapws.H2ONH3()


# fugacity in funciton of ptx
def fug(p, t, x):
    def hptx(rho):
        rho = max(1e-5, rho)
        y = mix._prop(rho, t, x)["P"]-p
        return y

    rhoo = (1-x)*iapws.IAPWS95._Liquid_Density(t) + x*NH3._Liquid_Density(t)
    rho = fsolve(hptx, rhoo)
    prop = mix._prop(rho, t, x)
    return prop["fugH2O"], prop["fugNH3"]


def _Bubble_T(P, zi):
    """Calculation Bubble Point Temperature"""
    # Initial estimation of temperature
    T = 0
    for xi, cmp in zip(zi, (iapws.IAPWS95, NH3)):
        Ti = cmp.Tc/(1-3*log(P/cmp.Pc)/(log(10)*(7+7*cmp.f_acent)))
        T += Ti*xi

    # Initial estimation of K using inverted Wilson correlation
    Ki = []
    for c in (iapws.IAPWS95, NH3):
        Ki.append(c.Pc/P*exp(5.37*(1.+c.f_acent)*(1.-c.Tc/T)))

    yi = [k*x for k, x in zip(Ki, zi)]
    ym = sum(yi)
    yi = [y/ym for y in yi]

    def f(zi, Ki):
        val = 0
        for z, ki in zip(zi, Ki):
            val += z*ki
        return val - 1

    c = 0
    while True:
        c += 1
        try:
            tital = fug(P, T, zi[1])
            titav = fug(P, T, yi[1])
        except OverflowError:
            T = None
            find = False
            break

        Ki = [l/v for l, v in zip(tital, titav)]
        if abs(f(zi, Ki)) <= 1e-5:
            find = 1
            break

        yi = [k*x for k, x in zip(Ki, zi)]
        ym = sum(yi)
        yi = [y/ym for y in yi]

        # Estimating next temperature value from [2]_, eq 11
        s = 0
        for cmp, x, k in zip((iapws.IAPWS95, NH3), zi, Ki):
            s += (1+cmp.f_acent)*cmp.Tc*x*k
        ed = T/5.373/s
        T *= (1-ed*f(zi, Ki))

        if c > 100:
            print("reach limit iteration count")
            find = 0
            break

    return T, find

P = 0.02
for x in arange(0, 1.1, 0.1):
    T, success = _Bubble_T(P, ((1-x), x))
    if not success:
        continue

    if x == 0:
        print("Pure water", iapws.IAPWS95(P=P, x=0).T-273.15)
        print(f"x={x}", T-273.15)
    elif x == 1:
        print(f"x={x}", T-273.15)
        print("Pure Ammonia", NH3(P=P, x=0).T-273.15)
    else:
        print(f"x={x}", T[0]-273.15)

Really I'm not sure the result are correct, differ something in the side values... Furthermore the convergence is really bad, and the algorithm is fairly inefficient, and this is the reason the equilibrium routine isn't implement.

Cheers, Juanjo

[bleves91]

thank you so much, I have something to start with. if I will get something I will share with you. I need this code to compute ammonia + water hybrid heat pump. they will be very popular in the next few years Il giorno lun 4 set 2023 alle ore 18:30 Juan José Gómez Romera < ***@***.***> ha scritto:

…

Hi, excuse me for my late response Try with this code from math import log, exp from numpy import arange from scipy.optimize import fsolve import iapws from iapws.ammonia import NH3 mix = iapws.H2ONH3() # fugacity in funciton of ptx def fug(p, t, x): def hptx(rho): rho = max(1e-5, rho) y = mix._prop(rho, t, x)["P"]-p return y rhoo = (1-x)*iapws.IAPWS95._Liquid_Density(t) + x*NH3._Liquid_Density(t) rho = fsolve(hptx, rhoo) prop = mix._prop(rho, t, x) return prop["fugH2O"], prop["fugNH3"] def _Bubble_T(P, zi): """Calculation Bubble Point Temperature""" # Initial estimation of temperature T = 0 for xi, cmp in zip(zi, (iapws.IAPWS95, NH3)): Ti = cmp.Tc/(1-3*log(P/cmp.Pc)/(log(10)*(7+7*cmp.f_acent))) T += Ti*xi # Initial estimation of K using inverted Wilson correlation Ki = [] for c in (iapws.IAPWS95, NH3): Ki.append(c.Pc/P*exp(5.37*(1.+c.f_acent)*(1.-c.Tc/T))) yi = [k*x for k, x in zip(Ki, zi)] ym = sum(yi) yi = [y/ym for y in yi] def f(zi, Ki): val = 0 for z, ki in zip(zi, Ki): val += z*ki return val - 1 c = 0 while True: c += 1 try: tital = fug(P, T, zi[1]) titav = fug(P, T, yi[1]) except OverflowError: T = None find = False break Ki = [l/v for l, v in zip(tital, titav)] if abs(f(zi, Ki)) <= 1e-5: find = 1 break yi = [k*x for k, x in zip(Ki, zi)] ym = sum(yi) yi = [y/ym for y in yi] # Estimating next temperature value from [2]_, eq 11 s = 0 for cmp, x, k in zip((iapws.IAPWS95, NH3), zi, Ki): s += (1+cmp.f_acent)*cmp.Tc*x*k ed = T/5.373/s T *= (1-ed*f(zi, Ki)) if c > 100: print("reach limit iteration count") find = 0 break return T, find P = 0.02 for x in arange(0, 1.1, 0.1): T, success = _Bubble_T(P, ((1-x), x)) if not success: continue if x == 0: print("Pure water", iapws.IAPWS95(P=P, x=0).T-273.15) print(f"x={x}", T-273.15) elif x == 1: print(f"x={x}", T-273.15) print("Pure Ammonia", NH3(P=P, x=0).T-273.15) else: print(f"x={x}", T[0]-273.15) Really I'm not sure the result are correct, differ something in the side values... Furthermore the convergence is really bad, and the algorithm is fairly inefficient, and this is the reason the equilibrium routine isn't implement. Cheers, Juanjo — Reply to this email directly, view it on GitHub <#68 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJVTACCCC6FEKCVU67A6RIDXYX62TANCNFSM6AAAAAAZWWC6OM> . You are receiving this because you authored the thread.Message ID: ***@***.***>