Are there any plans to implement Inside-out flash computing?

[ysyecust]

I have just come into contact with this project and would like to ask if flash computing currently implements inside-out mode calculations. If not, I can provide relevant help in the future.

[longemen3000]

I've heard about it, what are the main advantages over a michelsen's implementation? (successive substitution in fugacity + Gibbs minimization)

[ysyecust]

First of all, from the perspective of versatility, it is the most widely used flash calculation method in spen plus and Aspen hysys. Secondly, it has significant advantages in calculation speed and solution stability. Moreover, due to the characteristics of inner and outer calculation separation, many parallel optimizations Strategies can be added effectively; dealing with this type of problem at the same time can reduce the number of nested iterations as much as possible from an equation perspective. This is a relevant reference https://pubs.acs.org/doi/10.1021/acs.iecr.6b03956

[pw0908]

Interesting! This reminds me a little bit of the HELD algorithm, but only for two-phase region. Ill be curious to see how it compares to the Michelsen algorithm since that's usually the gold-standard.

I don't know if you want to give this a go @longemen3000, but @ysyecust, if you want to implement this, you can take a look at our implementation of the DETPFlash algorithm. If you have any trouble with the format, we'll be happy to help!

[longemen3000]

for T-P flash, it is something like this? it runs, but i suppose it is wrong (it does not converge to the real values of the flash model)

function inside_out_tp_flash(model, p, T, n, K0)
    z = n ./ sum(n)
    #0. input variables
    ε_out = 1e-10 #
    ε_in = 1e-10
    _1 = one(Base.promote_eltype(model,p,T,z,K0))
    
    #algorithm 4, 10.1021/acs.iecr.6b03956
    
    #1., 2. Solve 11) for α assuming Raoult's law
    α = _1*Clapeyron.rachfordrice(K0, z)
    
    #3. calculate x,y assuming Raoult's law
    x = Clapeyron.rr_flash_liquid!(fill(_1,length(z)),K0,z,α)
    y = Clapeyron.rr_flash_vapor!(fill(_1,length(z)),K0,z,α)
    x ./= sum(x)
    y ./= sum(y)
    #4. calculate k from eos
    lnϕx, volx = Clapeyron.lnϕ(model, p, T, x; phase=:l)
    lnϕy, voly = Clapeyron.lnϕ(model, p, T, y; phase=:v)

    k = exp.(lnϕx .- lnϕy)
    #5. set kb to 1
    kb = one(α)
    kb0 = 2*eps(kb)
    #6. calculate u (eq 14)
    u = log.(k ./ kb)
    û = similar(u)
    pᵢ = similar(u)
    kk = similar(u)
    #7. set R <- α
    R = α
    L = _1 - α
    V = α

    #8. set ΩPT <- 2*ε_out
    ΩPT = _1*2*ε_out
    ΨPT = _1*Inf
    #9 outer loop
    outer_count,inner_count = 0,0
    max_inner,max_outer = 100,100
    
    while ΩPT > ε_out
        outer_count += 1
        outer_count == max_outer && break
        inner_count = 0
        
        #10. inner loop
        while abs(ΨPT) > ε_in
            inner_count += 1
            inner_count == max_inner && break
            #11. calculate pᵢ,x,y,L from 22), 25) - 27)
            pᵢ .=  x ./ (1 .- R .+ kb0 .* R .* exp.(u))
            ∑pᵢ = sum(pᵢ)
            L = ∑pᵢ*(1 - R)
            x .= pᵢ
            x ./= ∑pᵢ
            V = 1 - L
            y .= exp.(u) .* pᵢ
            ∑euᵢpᵢ = sum(y)
            y  ./= sum(y)
            
            #12. calculate ΨPT from 53)
            ΨPT = mid(R, ∑pᵢ - kb*∑euᵢpᵢ, R - 1)
            kk .= y ./ x
            
            #13. assume new value for R (??)
            R = kb*V/(kb*V + kb0*L)
        end
       
        #15. calculate α from 6)
        #α = 1 - L

        #16. update K using EoS
        lnϕx, volx = Clapeyron.lnϕ!(lnϕx, model, p, T, x; phase=:l)
        lnϕy, voly = Clapeyron.lnϕ!(lnϕy, model, p, T, y; phase=:v)
        k .= exp.(lnϕx .- lnϕy)
    
        #17. calculate u (eq 14)
        û .= log.(k ./ kb)

        #18. calculate ΩPT from 54)
        ΩPT = Clapeyron.dnorm(u, û, Inf)
        #19. u <- û
        u .= û
    end
    return V,x,y
end

[pw0908]

Where did you get the mid() function?

[longemen3000]

sorry, this is mid: mid(a,b,c) = max(min(a,b),min(max(a,b),c))

[pw0908]

Made a few small changes (commented) but seems to be working now:

mid(a,b,c) =  max(min(a,b),min(max(a,b),c))

function inside_out_tp_flash(model, p, T, n, K0)
    z = n ./ sum(n)
    #0. input variables
    ε_out = 1e-10 #
    ε_in = 1e-10
    _1 = one(Base.promote_eltype(model,p,T,z,K0))
    
    #algorithm 4, 10.1021/acs.iecr.6b03956
    
    #1., 2. Solve 11) for α assuming Raoult's law
    α = _1*Clapeyron.rachfordrice(K0, z)
    
    #3. calculate x,y assuming Raoult's law
    x = Clapeyron.rr_flash_liquid!(fill(_1,length(z)),K0,z,α)
    y = Clapeyron.rr_flash_vapor!(fill(_1,length(z)),K0,z,α)
    x ./= sum(x)
    y ./= sum(y)
    #4. calculate k from eos
    lnϕx, volx = Clapeyron.lnϕ(model, p, T, x; phase=:l)
    lnϕy, voly = Clapeyron.lnϕ(model, p, T, y; phase=:v)

    k = exp.(lnϕx .- lnϕy)
    #5. set kb to 1
    kb = one(α)
    kb0 = 2*eps(kb)
    #6. calculate u (eq 14)
    u = log.(k ./ kb)
    û = similar(u)
    pᵢ = similar(u)
    kk = similar(u)
    #7. set R <- α
    R = α
    L = _1 - α
    V = α

    #8. set ΩPT <- 2*ε_out
    ΩPT = _1*2*ε_out
    ΨPT = _1*Inf
    #9 outer loop
    outer_count,inner_count = 0,0
    max_inner,max_outer = 100,100
    
    while ΩPT > ε_out
        outer_count += 1
        outer_count == max_outer && break
        inner_count = 0
        
        #10. inner loop
        while abs(ΨPT) > ε_in
            inner_count += 1
            inner_count == max_inner && break
            #11. calculate pᵢ,x,y,L from 22), 25) - 27)
            pᵢ .=  z ./ (1 .- R .+ kb .* R .* exp.(u)) # Change
            ∑pᵢ = sum(pᵢ)
            L = ∑pᵢ*(1 - R)
            x .= pᵢ
            x ./= ∑pᵢ
            V = 1 - L
            y .= exp.(u) .* pᵢ
            ∑euᵢpᵢ = sum(y)
            y  ./= sum(y)
            
            #12. calculate ΨPT from 53)
            ΨPT = mid(R, ∑pᵢ - kb*∑euᵢpᵢ, R - 1)
            kk .= y ./ x
            
            #13. assume new value for R (??)
            R = kb*V/(kb*V + L) # Changed
        end
       
        #15. calculate α from 6)
        # α = 1 - L

        #16. update K using EoS
        lnϕx, volx = Clapeyron.lnϕ!(lnϕx, model, p, T, x; phase=:l)
        lnϕy, voly = Clapeyron.lnϕ!(lnϕy, model, p, T, y; phase=:v)
        k .= exp.(lnϕx .- lnϕy)
    
        #17. calculate u (eq 14)
        û .= log.(k ./ kb)

        #18. calculate ΩPT from 54)
        ΩPT = Clapeyron.dnorm(u, û, Inf)
        #19. u <- û
        u .= û
    end
    return V,x,y
end