分离工程泡露点计算

最新推荐文章于 2024-07-25 18:55:58 发布

置顶 MyZinG996

最新推荐文章于 2024-07-25 18:55:58 发布

阅读量926

点赞数 4

文章标签： python

本文链接：https://blog.csdn.net/Taylorzing/article/details/129889690

版权

A mixture with the composition given below is fed to a rectifying column at the operating pressure of 0.1MPa. Calculate the feeding temperature at bubble | dew point.

Component	butane	pentane	n-hexane	heptane	octane	total
Mole fraction	0.05	0.17	0.65	0.10	0.03	1

Antoine equation:

p:bar( $10^{5}$ Pa）; T:K

$K_{i}=\frac{p^{0}}{p^{s}}$

$P^s$ :混合总压 :纯组分的饱和蒸汽压

输入基本参数：

import numpy as np
from sympy import *
from scipy.misc import derivative
P = 0.1  # 0.1Mpa=0.1bar
x = [0.05, 0.17, 0.65, 0.10, 0.03]
a = [9.0580, 9.2131, 9.2164, 9.2535, 9.3224]
b = [2154.90, 2477.07, 2697.55, 2911.32, 3120.29]
c = [-34.42, -39.94, -48.8, -56.51, -63.63]
t = symbols("t")

$\sum y _{i}=\sum K_{_{i}}\cdot x_{_{i}}$

ysum = 0
for i in range(len(x)):
    # f = exp(a[i] + b[i] / t + c[i] * P) * x[i]
    f = (exp(a[i] - b[i] / (c[i] + t)) / P) * x[i]
    print(f)
    ysum += f
print(ysum)

迭代方法：

1.牛顿法

def Newton(x0, f, epsilon=1 / 10 ** 5, iternum=100):  # 初值，精度要求，最大迭代次数,迭代函数
    xk_1 = x0
    fdx = diff(f, t)
    for i in range(iternum):
        # fx = f(xk_1)
        # fdx = derivative(f, xk_1, dx=1e-6)
        fx = f.subs(t, xk_1)

        fdx = fdx.subs(t, xk_1)
        if fdx != 0:
            xk = xk_1 - fx / fdx
            print("第", i + 1, "次迭代 ", "xk=", xk, "  xk-1=", xk_1, "  |xk - xk-1|=", abs(xk - xk_1));
            if abs(xk - xk_1) < epsilon:
                return xk
            else:
                xk_1 = xk
        else:
            break
    print("方法失败")
    return 0

2.弦截法

def xqf(f, a=300, b=499, epsilion=1 / 10 ** 4, iter=3000):

    xa, xb = a, b
    for i in range(iter):
        fa = f.subs(t, xa)
        fb=f.subs(t, xb)
        # fb = f.evalf(f.subs(t, xb))
        x_1 = xb - (fb * (xb - xa) / (fb - fa))
        print(fa,fb)
        print("第", i + 1, "次迭代 ", "xa=", xa, "  x_1=", x_1, "  |xa - xb|=", abs(xa - xb))
        if abs(xb - xa) < epsilion:
                break
        else:
            xa, xb = xb, x_1
    return x_1

3 .Richmond法

def richmond(f, x0=100, epsilon=1 / 10 ** 5, iternum=50):
    global count
    xk_1 = x0
    fdx = diff(f, t)
    fdx2 = diff(f, t, 2)
    xmin = []
    print("该函数第%d次递归" % count)
    count += 1
    for i in range(iternum):

        fx = f.subs(t, xk_1)
        fdx = fdx.subs(t, xk_1)
        fdx2 = fdx2.subs(t, xk_1)

        if fdx != 0:
            xk = xk_1 - 2 / ((2 * fdx / fx) - (2 * fdx2 / fdx))
            print("第", i + 1, "次迭代 ", "xk=", xk, "  xk-1=", xk_1, "  |xk - xk-1|=", abs(xk - xk_1))
            if abs(xk - xk_1) < 2:
                xmin.append(xk_1)
            if abs(xk - xk_1) < epsilon:
                return xk
            else:
                xk_1 = xk
        else:
            break
        # print(xmin)
    xnew = min(xmin)  # 返回最小值
    print(xnew)
    return xnew

个人总结：

对于迭代方法初值x0的取值陷入不收敛或局部最优，可以预估其初值（二元函数较为容易————我采用画图）；建议限定条件避免此情况发生（没学过，感觉可以这么做）。
泡露点计算均为完全理想体系（活度系数，逸度系数难求，还要内层迭代试差麻烦，能力有限）

参考化工热力学：完全非理想系的求法

$K_{i}=\frac{p_{i}^{0}\phi _{0}^{i}\gamma _{i}}{p\phi _{i}^{V}}$

总计算框图

Aspen获得数据最后三行为：Tc：临界温度 Pc：临界压力 偏心因子：ω（OMEGA）

活度系数求法（Wilson（公式去问aspen））

$ln\lambda _{i}=1-ln(\sum_{j}^{}A_{ij}x_j)-\sum_{j}\frac{A_{ji}x_{j}}{\sum_{k}A_{jk}x_k}$

逸度系数

$ln\phi_{i}=\frac{P_{r}}{T_{r}}(B^{0}+B^{1})$

$T_{r}=\frac{{T}}{T_{c}}$ $P_{r}=\frac{P}{P_{r}}$
def fugacity_coefficient():
    tc = [425.12, 269.7, 507.6, 540.2]  # K
    pc = [37.96, 33.7, 30.25, 27.4]  # bar(10^5)
    w = [0.200164, 0.251206, 0.301261, 0.349469]  # 偏心因子ω
    pr, tr, fc, B0, B1 = [], [], [], [], []
    for i in range(len(y) - 1):
        pr.append(P / pc[i])
        tr.append(t / tc[i])
        B0.append(0.083 - 0.422 / pow(tr[i], 1.6))
        B1.append(0.139 - 0.172 / pow(tr[i], 4.2))
        fc.append(sy.exp(pr[i] / tr[i])*(B0[i] + w[i] * B1[i]))
    return fc

可以尝试其他题目：

其他例题：

A mixture contains 5 mol% ethane, 30 mol% propane, and 65 mol% n-butane. The equilibrium constants of these components under the operating pressure can be calculated by the following equations. Calculate the bubble point temperature.

ethane: K=0.13333t+4.6667 ; propane:K=0.6667t+1.13333;

n-butane:K=0.02857t+0.08571 (t, ℃)

The composition and equilibrium constants of a vapor mixture are as follows:

Component

A

B

C

Mole fraction

0.35

0.2

0.45

Ki (T，℃；p，MPa）
0.15T/P
0.02T/P

0.01T/P

Calculate the dew point temperature at the pressure of 2MPa (the error criterion can be 0.001)

A 100kmol/h feed consisting of 10, 20, 30 and 40 mol% of propane (1), n-butane (2) n-pentane (3), and n-hexane (4), respectively, enters a distillation column at 366.5K and 689.5kPa. Assuming equilibrium, calculate mole fraction of the feed enters as vapor, and the liquid and vapor compositions.

（此题最难，没学懂且程序复杂，希望有大佬解出【求代码，评论区见】）

A mixture consists of 40 mol%, 30 mol% and 30 mol% of n-butane, n-pentane and n-hexane, respectively. The total pressure is 1.013×103 kPa.

(1) Calculate the bubble point of the mixture and the vapor composition assuming equilibrium.

(2) Calculate mole fraction of the feed enters as vapor, and the liquid and vapor compositions at 122 °C and 1.013 × 103 kPa.

x = [0.05, 0.3, 0.65]
t = symbols("t")
K = [0.1333 * t + 4.6667, 0.6667 * t + 1.3333, 0.02857 * t + 0.08571]
Y = -1
for i in range(0, len(x)):
    y = x[i] * K[i]
    Y += y
def Newton(x0, epsilon, f,iternum=100):  # 初值，精度要求，最大迭代次数,迭代函数
    xk_1 = x0
    fdx = diff(f, t)
    for i in range(iternum):
        # fx = f(xk_1)
        # fdx = derivative(f, xk_1, dx=1e-6)
        fx=f.subs(t,xk_1)
        fdx=fdx.subs(t,xk_1)
        if fdx != 0:
            xk = xk_1 - fx / fdx
            print("第", i + 1, "次迭代 ", "xk=", xk, "  xk-1=", xk_1, "  |xk - xk-1|=", abs(xk - xk_1));
            if abs(xk - xk_1) < epsilon:
                return xk
            else:
                xk_1 = xk
        else:
            break
    print("方法失败")
    return 0
Newton(30,1/10**4,f=Y)

import matplotlib.pyplot as plt
import numpy as np
from numpy import linspace
from sympy import *
from scipy.misc import derivative

P = 1  # 0.1Mpa=1bar
y = [0.05, 0.17, 0.65, 0.10, 0.03]
# x = [0.05, 0.17, 0.65, 0.10, 0.03]
a = [9.0580, 9.2131, 9.2164, 9.2535, 9.3224]
b = [2154.90, 2477.07, 2697.55, 2911.32, 3120.29]
c = [-34.42, -39.94, -48.8, -56.51, -63.63]
t = symbols("t")


def Newton(x0, f, epsilon=1 / 10 ** 5, iternum=100):  # 初值，精度要求，最大迭代次数,迭代函数
    xk_1 = x0
    fdx = diff(f, t)
    for i in range(iternum):
        # fx = f(xk_1)
        # fdx = derivative(f, xk_1, dx=1e-6)
        fx = f.subs(t, xk_1)

        fdx = fdx.subs(t, xk_1)
        if fdx != 0:
            xk = xk_1 - fx / fdx
            print("第", i + 1, "次迭代 ", "xk=", xk, "  xk-1=", xk_1, "  |xk - xk-1|=", abs(xk - xk_1));
            if abs(xk - xk_1) < epsilon:
                return xk
            else:
                xk_1 = xk
        else:
            break
    print("方法失败")
    return 0


xsum = 0
for i in range(len(y)):
    # f = exp(a[i] + b[i] / t + c[i] * P) * x[i]
    f = y[i]/(exp(a[i] - b[i] / (c[i] + t)) / P)
    print(f)
    xsum += f

print(xsum)


def dawfunc():
    tt = linspace(300, 450)
    yy = []
    for t in tt:
        yy.append(xsum.subs(symbols("t"), t))
        # func.evalf(subs={'x2': 6})
        # yy.append(ysum.evalf(subs={'t': t}))

    plt.plot([300,450], [1,1])
    plt.plot(tt, yy)
    plt.show()


# Newton(x0=1000,f=ysum-1,epsilon=1 / 10 ** 5, iternum=1000)
def xqf(f, a=300, b=499, epsilion=1 / 10 ** 4, iter=3000):
    xa, xb = a, b
    for i in range(iter):
        fa = f.subs(t, xa)
        fb=f.subs(t, xb)
        # fb = f.evalf(f.subs(t, xb))
        x_1 = xb - (fb * (xb - xa) / (fb - fa))
        print(fa,fb)
        print("第", i + 1, "次迭代 ", "xa=", xa, "  x_1=", x_1, "  |xa - xb|=", abs(xa - xb))
        if abs(xb - xa) < epsilion:
                break
        else:
            xa, xb = xb, x_1
    return x_1


xqf(f=xsum - 1, a=330, b=411, epsilion=1 / 10 ** 10, iter=300)
Newton(x0=320,f=xsum-1,epsilon=1 / 10 ** 5, iternum=1000)
dawfunc()

引用文献：一种新的泡点计算方法https://kns.cnki.net/kcms2/article/abstract?v=3uoqIhG8C44YLTlOAiTRKjkpgKvIT9NkGsvn6cq9Bo2RsFREO0F3IAc2j6lsyDNGl5xjeV9fG5aEGsOvsz2TEbiNonMKNexZ&uniplatform=NZKPT

MyZinG996

关注

4
点赞
踩
4

收藏

觉得还不错? 一键收藏
3
评论
分离工程泡露点计算

对于迭代方法初值x0的取值陷入不收敛或局部最优，可以预估其初值（二元函数较为容易————我采用画图）；建议限定条件避免此情况发生（没学过，感觉可以这么做）。泡露点计算均为完全理想体系（活度系数，逸度系数难求，还要内层迭代试差麻烦，能力有限）参考化工热力学。
复制链接

扫一扫