利用python图解法求理论板层数

参考文章链接

理论板层数的计算是精馏塔设计型计算的主要内容，通过塔板效率便可确定塔的有效高度。计算理论板层数的基本步骤是：规定分离要求，确定操作条件，利用平衡关系和操作关系计算所需理论板层数。工程确定理论板层数可采取逐板及算法，图解法和简捷法，本文将用图解法来计算理论板层数。

  • 分离要求：原料液流量 100kmol/h 组成为0.44 镏出液组成 0.975釜残成 0.0235
  • 操作条件：露点进料 相对挥发度 2.47

1 导入所需模块

1
2
3
4
5
import numpy as np
import sympy
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

2 根据分离要求和操作条件确定参数

1
2
3
4
5
arguement = [0.44,0.975,0.0235, None,0]
#挥发度
alpha = 2.47
F = 100
xF,xD,xW,R,q = arguement

3 分别定义 气液平衡、操作线、q线方程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 平衡线方程
#用x表示y的平衡线方程
def yp(x):
    y = alpha*x / (1+(alpha-1)*x)
    return y
#用y表示x的平衡线方程
def xp(y):
    x = y/(alpha - (alpha -1)*y)
    return x
# 精馏段操作方程
def yj(x):
    y = R/(R+1)*x + xD/(R+1)
    return y
# q线方程
def yq(x, q):
    if q == 1:
        pass
    else:
        return q/(q-1)*x -1/(q-1)*xF

此处先利用线代的方法求出D,W再利用提馏段操作方程

1
2
3
4
5
6
7
8
A = sympy.Matrix([[1,1],[xD,xW]])
b = sympy.Matrix([F,xF*F])
D,W = A.solve(b)
L = R*D
#   提馏段操作方程
def yt(x):
    y = (L + q*F)/(L + q*F -W)*x - W/(L + q*F -W)*xW
    return y

计算回流比
1
2
3
4
R_min = 1/(alpha - 1)*(alpha*xD/xF - (1-xD)/1-xF) - 1
R = R_min * 1.5
#当为饱和液体进料时
R_min = 1/(alpha - 1)*(alpha/xF - alpha*(1-xD)/1-xF)

4 逐板计算，求解每个塔板的平衡情况

1
2
3
4
5
6
7
8
9
10
11
yn = [xD]
xn = []

while xp(yn[-1])>xW:
    xn.append(xp(yn[-1]))
    if xn[-1]>xq:
        yn.append(yj(xn[-1]))
    else:
        yn.append(yt(xn[-1]))
else:
    xn.append(xp(yn[-1]))

5 图解法计算理论塔板数的图示数据

1
2
3
4
5
6
7
8
9
10
xnt = [xD]
ynt = [xD]
for n,i in enumerate(xn):
        xnt.append(i)
        ynt.append(yn[n])
        xnt.append(i)
        if i >= xq:
            ynt.append(yj(i))
        else:
            ynt.append(yt(i))

6 作图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ax.plot(x,x,ls=':',label='对角线') #画线
ax.plot(x,yq1,label='q线')
ax.plot(x,yp1,label='平衡线')
ax.plot(x,yj1,label='精馏操作线')
ax.plot(x,yt1,label='提馏操作线')
ax.plot(xn,yn,label='塔板操作平衡点',ls=':',marker='+',markersize=10)
ax.plot(xnt,ynt,label='图解法—理论塔板',ls=':')

ax.plot(xD,xD,marker='.',markersize=10) #画点
ax.plot(xW,xW,marker='.',markersize=10)
ax.plot(xq,yq,marker='.',markersize=10)

ax.annotate('W点',xy=(xW,xW),xytext=(xW+0.06,xW),arrowprops=dict(arrowstyle='->'))#画点的注释
ax.annotate('D点',xy=(xD,xD),xytext=(xD,xD-0.06),arrowprops=dict(arrowstyle='->'))
ax.annotate("Q点", xy=(xq, yq), xytext=(xq, yq - 0.05), arrowprops=dict(arrowstyle="->"))

如图所示

7 完整代码如下所示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#导入模块
import numpy as np
import sympy
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
# 参数e条件:原料液流量 100kmol/h  组成为0.44   镏出液组成 0.975
# 釜残液组成 0.0235  露点进料  相对挥发度 2.47
arguement = [0.44,0.975,0.0235, None,0]
#挥发度
alpha = 2.47
F = 100
xF,xD,xW,R,q = arguement

# q线方程
def yq(x, q):
    if q == 1:
        pass
    else:
        return q/(q-1)*x -1/(q-1)*xF
# 平衡线方程
#用x表示y的平衡线方程
def yp(x):
    y = alpha*x / (1+(alpha-1)*x)
    return y
#用y表示x的平衡线方程
def xp(y):
    x = y/(alpha - (alpha -1)*y)
    return x
#回流比的计算
# 饱和蒸汽进料---q==0
R_min = 1/(alpha - 1)*(alpha*xD/xF - (1-xD)/1-xF) - 1
R = R_min * 1.5
#当为饱和液体进料时
R_min = 1/(alpha - 1)*(alpha/xF - alpha*(1-xD)/1-xF)
parameter = [xF,xD,xW,R,q]
print('R:{}'.format(R))

# 精馏段操作方程
def yj(x):
    y = R/(R+1)*x + xD/(R+1)
    return y

# 提馏段操作方程
# 需要先求解D和W,为此，使用sympy求解
# 使用线性代数的方法连解 F=D+W 和  F*xF=DxD+WxW
# A = sympy.Matrix([[1,1],[D,W]]) 创建了一个 2x2 的矩阵
# b = sympy.Matrix([F,F*xf]) 创建了一个列向量 b，
A = sympy.Matrix([[1,1],[xD,xW]])
b = sympy.Matrix([F,xF*F])
D,W = A.solve(b)
L = R*D
#   提馏段操作方程
def yt(x):
    y = (L + q*F)/(L + q*F -W)*x - W/(L + q*F -W)*xW
    return y

# 准备数据，q线，平衡线，精馏操作线，提馏操作线
x = np.linspace(0,1,50)
yq1 = [yq(xi, q) for xi in x]
yp1 = yp(x)
yj1 = yj(x)
yt1 = yt(x)

# 确定Q点
xq = ((R+1)*xF + (q-1)*xD)/(R+q)
yq = (xF*R + q*xD)/(R+q)

# 逐板计算，求解每个塔板的平衡情况
yn = [xD]
xn = []

while xp(yn[-1])>xW:
    xn.append(xp(yn[-1]))
    if xn[-1]>xq:
        yn.append(yj(xn[-1]))
    else:
        yn.append(yt(xn[-1]))
else:
    xn.append(xp(yn[-1]))
    print('N_T = {}'.format(len(xn)))
    count = sum(1 for x in xn if x > xq)
    print(f"大于 {xq} 的元素个数是：{count}")


# 输出塔板上的平衡点，经过四舍五入的结果
yn_r = [round(i,3) for i in yn]
xn_r = [round(i,3) for i in xn]

print('塔板上的平衡点，经过四舍五入的结果:')
print('yn_r={}'.format(yn_r))
print('xn_r={}'.format(xn_r))

# 图解法计算理论塔板数的图示数据
xnt = [xD]
ynt = [xD]
for n,i in enumerate(xn):
        xnt.append(i)
        ynt.append(yn[n])
        xnt.append(i)
        if i >= xq:
            ynt.append(yj(i))
        else:
            ynt.append(yt(i))
# 画图
#基础设置
fig, ax = plt.subplots(1, 1, figsize=(9, 9))
ax.set_xlim(0,1)
ax.set_ylim(0,1)

#   画图
ax.plot(x,x,ls=':',label='对角线',color="black") #画线
ax.plot(x,yq1,label='q线')
ax.plot(x,yp1,label='平衡线',color="black")
ax.plot(x,yj1,label='精馏操作线')
ax.plot(x,yt1,label='提馏操作线')
ax.plot(xn,yn,label='塔板操作平衡点',ls=':',marker='+',markersize=10)
ax.plot(xnt,ynt,label='图解法—理论塔板',ls=':',color="red")

ax.plot(xD,xD,marker='.',markersize=10) #画点
ax.plot(xW,xW,marker='.',markersize=10)
ax.plot(xq,yq,marker='.',markersize=10)

ax.annotate('W点',xy=(xW,xW),xytext=(xW+0.06,xW),arrowprops=dict(arrowstyle='->'))#画点的注释
ax.annotate('D点',xy=(xD,xD),xytext=(xD,xD-0.06),arrowprops=dict(arrowstyle='->'))
ax.annotate("Q点", xy=(xq, yq), xytext=(xq, yq - 0.05), arrowprops=dict(arrowstyle="->"))
# 获取坐标轴
# ax = plt.gca()
# 设置坐标轴顶部线条粗细
ax.spines["top"].set_linewidth(2)
# 设置坐标轴底部线条粗细
ax.spines["bottom"].set_linewidth(2)
# 设置坐标轴左侧线条粗细
ax.spines["left"].set_linewidth(2)
# 设置坐标轴右侧线条粗细
ax.spines["right"].set_linewidth(2)
# ax.grid()
# 图中显示所需理论板数
ax.text(x=0.6, y=0.4, s="所需理论板数：%d" % (len(xn) - 1))
ax.text(x=0.6, y=0.3, s="进料板是第%d层" % (count+1))
ax.legend()
plt.show()