多效蒸发在化工、轻工、食品工业及海水淡化等行业中得到了广泛应用。许多研究者对常规并流多效蒸发系统模拟计算进行了研究, 大多数求解是采用试差法[1 ] 、牛顿迭代法[2 ] 和Broyden法[3 ] 。众所周知, 蒸发大量水分时, 要消耗大量的加热蒸汽。为了减少加热蒸汽消耗, 提高生蒸汽的经济程度, 除采用多效蒸发以外, 还可采用冷凝水闪蒸、引出额外蒸汽预热原料液等措施。然而, 带有冷凝水闪蒸和引出额外蒸汽的复杂并流多效蒸发系统的数学模型及其求解方法未见报道。笔者首次建立了复杂并流多效蒸发系统的数学模型, 提出一种新算法——迭代法结合矩阵法求解模型并给出算法框图, 用Visual Basic 5.0语言编制了该算法的计算程序。
《1 数学模型》
1 数学模型
《1.1 工艺流程》
1.1 工艺流程
有冷凝水闪蒸和额外蒸汽引出的复杂并流多效蒸发系统的流程如图1所示。为了进一步利用冷凝水的热能, 图1的流程中将从前一效冷凝水闪蒸器流出的温度较高的冷凝水引到后一效温度较低的冷凝水闪蒸器中继续闪蒸。另外, 设除末效外, 各效均有冷凝水闪蒸和引出额外蒸汽。系统总效数为n , 图1中只示意画到第i 效。
《1.2 系统物料衡算》
1.2 系统物料衡算
假定溶质不挥发, 对任一效 (第i 效) 的溶质进行衡算, 有
F 0 w 0 = F i w i = ( F 0 − W 1 − ⋯ − W i ) w i ( 1 ) F 0 w 0 = F i w i = ( F 0 - W 1 - ⋯ - W i ) w i ( 1 )
W = ∑ n i = 1 W i = F 0 ( 1 − w 0 / w n ) ( 2 ) W = ∑ n i = 1 W i = F 0 ( 1 - w 0 / w n ) ( 2 )
《1.3 系统热量衡算》
1.3 系统热量衡算
考虑热损失及浓缩热, 对图1所示的复杂并流多效蒸发系统第i 效先按无冷凝水闪蒸和额外蒸汽引出的情况进行热量衡算可得[1 ]
W i = { α i D i + [ F 0 c 0 − c ∗ ( W 1 + ⋯ + W i − 1 ) ] β i } η i ( 3 ) W i = { α i D i + [ F 0 c 0 - c * ( W 1 + ⋯ + W i - 1 ) ] β i } η i ( 3 )
式中:η i 为第i 效的热利用系数, 对于一般溶液的蒸发, 可取η i =0.96~0.98, 对于浓缩热较大的物料, 可取η i =0.98-0.700 (w i -w i -1 ) [1 ] ;α i 、β i 分别为第i 效的蒸发系数和自蒸发系数, 其定义式如下
α i = ( H i − 1 − c ∗ t ′ i − 1 ) / ( H i − c ∗ t i ) β i = ( t i − 1 − t i ) / ( H i − c ∗ t i ) α i = ( Η i - 1 - c * t ′ i - 1 ) / ( Η i - c * t i ) β i = ( t i - 1 - t i ) / ( Η i - c * t i )
《图1》
图1 复杂并流多效蒸发系统流程示意图
Fig.1 The flow chart of complex cocurrent multi-effect evaporation
《图2》
现在对图1所示的复杂并流流程, 考虑有冷凝水闪蒸与额外蒸汽引出, 有
i = 1 时 ‚ D i = D 1 i ≥ 2 时 ‚ D i = W i − 1 − E i − 1 + G i − 1 } ( 4 ) i = 1 时 ‚ D i = D 1 i ≥ 2 时 ‚ D i = W i - 1 - E i - 1 + G i - 1 } ( 4 )
当i≥3时, 对第i-1效冷凝水闪蒸器进行热量衡算, 可得
⎛ ⎝ ⎜ ⎜ ∑ i − 1 i = 1 D i − ∑ i − 2 i = 1 G i ⎞ ⎠ ⎟ ⎟ c ∗ t ′ i − 2 = G i − 1 H i − 1 + ⎡ ⎣ ⎢ ⎢ ∑ i − 1 i = 1 ( D i − G i ) ⎤ ⎦ ⎥ ⎥ c ∗ t ′ i − 1 = G i − 1 ( H i − 1 − c ∗ t ′ i − 1 ) + ⎛ ⎝ ⎜ ⎜ ∑ i − 1 i = 1 D i − ∑ i − 2 i = 1 G i ⎞ ⎠ ⎟ ⎟ c ∗ t ′ i − 1 ( ∑ i - 1 i = 1 D i - ∑ i - 2 i = 1 G i ) c * t ′ i - 2 = G i - 1 Η i - 1 + [ ∑ i - 1 i = 1 ( D i - G i ) ] c * t ′ i - 1 = G i - 1 ( Η i - 1 - c * t ′ i - 1 ) + ( ∑ i - 1 i = 1 D i - ∑ i - 2 i = 1 G i ) c * t ′ i - 1
G i − 1 = ⎛ ⎝ ⎜ ⎜ ∑ i − 1 i = 1 D i − ∑ i − 2 i = 1 G i ⎞ ⎠ ⎟ ⎟ c ∗ ω i − 1 = ⎛ ⎝ ⎜ ⎜ D 1 + ∑ i − 1 i = 2 D i − ∑ i − 2 i = 1 G i ⎞ ⎠ ⎟ ⎟ c ∗ ω i − 1 ( 5 ) G i - 1 = ( ∑ i - 1 i = 1 D i - ∑ i - 2 i = 1 G i ) c * ω i - 1 = ( D 1 + ∑ i - 1 i = 2 D i - ∑ i - 2 i = 1 G i ) c * ω i - 1 ( 5 )
式中:ωi-1 为第i-1效冷凝水闪蒸器的闪蒸系数, 其定义式如下
ω i − 1 = ( t ′ i − 2 − t ′ i − 1 ) / ( H i − 1 + c ∗ t ′ i − 1 ) ω i - 1 = ( t ′ i - 2 - t ′ i - 1 ) / ( Η i - 1 + c * t ′ i - 1 )
利用式 (4) 的结果 (Di -Gi-1 =Wi-1 -Ei-1 ) 则式 (5) 可写成
G i − 1 = ( D 1 + W 1 + ⋯ + W i − 2 − E 1 − ⋯ − E i − 2 ) c ∗ ω i − 1 ( 6 ) G i - 1 = ( D 1 + W 1 + ⋯ + W i - 2 - E 1 - ⋯ - E i - 2 ) c * ω i - 1 ( 6 )
同理, 当i=2时, 对第1效冷凝水闪蒸器进行热量衡算, 有
G 1 = D 1 c ∗ ( t ′ 0 − t ′ 1 ) / ( H 1 − c ∗ t ′ 1 ) = D 1 c ∗ ω 1 ( 7 ) G 1 = D 1 c * ( t ′ 0 - t ′ 1 ) / ( Η 1 - c * t ′ 1 ) = D 1 c * ω 1 ( 7 )
《1.4 各效传热面积的计算》
1.4 各效传热面积的计算
A i = Q i / ( K i Δ t i ) = D i r i / ( K i Δ t i ) ( 8 ) A i = Q i / ( Κ i Δ t i ) = D i r i / ( Κ i Δ t i ) ( 8 )
式中:ri 为第i效加热蒸汽 (即i-1效二次饱和蒸汽) 的汽化潜热, 对第1效而言即为生蒸汽汽化潜热;Δ ti 为第i效有效传热温差, 当i=1时, Δ t1 =t′0 -t1 , 当i≥2时, Δ ti 应扣去Δ i 即Δ ti =t′i-1 -Δ i -ti , Δ i 值根据经验可取1℃[1 ] ;Di 根据式 (4) 中定义处理;Δ ti 按等面积原则[1 ] 进行分配
Δ t i = Δ t ( D i r i / K i ) / ⎛ ⎝ ⎜ ∑ n i = 1 D i r i / K i ⎞ ⎠ ⎟ Δ t = ∑ n i = 1 Δ t i = t ′ 0 − t ′ K − ∑ n i = 1 ( Δ ′ i + Δ ′′ i + Δ ′′′ i ) ⎫ ⎭ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ( 9 ) Δ t i = Δ t ( D i r i / Κ i ) / ( ∑ n i = 1 D i r i / Κ i ) Δ t = ∑ n i = 1 Δ t i = t ′ 0 - t ′ Κ - ∑ n i = 1 ( Δ ′ i + Δ ″ i + Δ ‴ i ) } ( 9 )
计算αi 、βi 、Δ ti 时均需各效溶液的沸点ti , 其值可按相平衡方程[1 ] (汽液温度关系式) 计算
t i = t ′ i + Δ ′ i + Δ ′′ i t i = t ′ i + Δ ′ i + Δ ″ i
式中:Δ ′i 用杜林方程[1 ] 或专用的关联式[4 ] 计算;Δ ″i 与各效蒸发压力Pi 及溶液液面高度有关[1 ] , Pi 可通过安托因方程与t′i 关联后计算。
αi \, βi 和ωi-1 的定义式均与饱和蒸汽的焓值Hi 有关。为编程计算方便, 将文献[1 ] 饱和水蒸汽表中的汽化潜热和焓与温度的关系回归成如下:
H i = 2 4 7 4 7 7 1 . 0 + 2 4 1 0 . 2 t ′ i − 3 . 8 t ′ 2 i r i = 2 4 6 6 9 0 4 . 9 − 1 5 8 4 . 3 t ′ i − 1 − 4 . 9 t ′ 2 i − 1 Η i = 2 4 7 4 7 7 1 . 0 + 2 4 1 0 . 2 t ′ i - 3 . 8 t ′ i 2 r i = 2 4 6 6 9 0 4 . 9 - 1 5 8 4 . 3 t ′ i - 1 - 4 . 9 t ′ i - 1 2
《2 模型求解》
2 模型求解
多效蒸发计算就是联立求解各效物料衡算式、热量衡算式、相平衡方程式及传热速率方程式等, 从而求出Di 、Wi 、wi 和Ai 等。上述算式实质是一个非线性方程组, 当效数比较多、方程组庞大时, 求解较困难。多数文献采用试差法[1 ] 、牛顿迭代法[2 ] 和Broyden 法[3 ] 求解。上述算法中Broyden 法效率较高, 但仍需取较多的初值, 需要较多次的迭代才会收敛。文献[3 ] 报道了用Broyden 法求解常规三效并流蒸发过程, 共12个方程, 需取12个初值, 经过19次迭代收敛。其中有些初值如D1 、Ai 、ti 、t′i 等凭经验选取较困难, 势必会影响到收敛的稳定性和迭代次数。本文建立的是有冷凝水闪蒸和额外蒸汽引出的复杂并流多效蒸发系统, 变量和方程的数目更多, 求解也更困难, 因而有必要寻求更有效的求解方法。
《2.1 系统物料衡算和热量衡算方程组的矩阵形式》
2.1 系统物料衡算和热量衡算方程组的矩阵形式
把式 (4) 代入式 (3) 可得各效蒸发水分量的计算式, 由式 (6) 和式 (7) 可得各效冷凝水闪蒸汽量的计算式, 由式 (2) 可得总蒸发水分量计算式, 上述各式可写成矩阵方程, 待求的未知量有2n个 (包括D1 , W1 , W2 , …, Wn , G1 , G2 , …, Gn-1 ) 。该矩阵是稀疏矩阵, 将其写成分块矩阵会使原矩阵结构简单运算更容易。分析该矩阵的结构, 将其分成八块处理比较适宜, 其具体结构如下:
[ A 1 A 2 A 3 A 4 ] × [ B 1 B 2 ] = [ C 1 C 2 ] [ A 1 A 3 A 2 A 4 ] × [ B 1 B 2 ] = [ C 1 C 2 ]
A 1 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ α 1 η 1 0 0 ⋮ 0 ⋮ 0 0 − 1 α 2 η 2 − c ∗ β 2 η 2 − c ∗ β 3 η 3 ⋮ − c ∗ β i η i ⋮ − c ∗ β n η n 1 − 1 α 3 η 3 − c ∗ β 3 η 3 ⋮ ⋯ ⋮ ⋯ ⋯ − 1 ⋮ − c ∗ β i η i α i η i ⋮ ⋯ 1 ⋱ − c ∗ β i η i ⋮ − c ∗ β n η n ⋯ − 1 ⋮ − c ∗ β n η n α n η n 1 ⋱ c ∗ β n η n 1 − 1 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ , A 2 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ c ∗ ω 1 c ∗ ω 2 ⋮ c ∗ ω i ⋮ c ∗ ω n − 1 c ∗ ω 2 ⋯ ⋯ ⋱ ⋯ ⋯ c ∗ ω i ⋯ ⋯ c ∗ ω n − 1 0 0 0 0 0 0 0 0 0 0 0 0 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ‚ A 3 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0 α 2 η 2 α 3 η 3 ⋱ α i η i ⋱ α n η n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ‚ A 4 = ⎡ ⎣ ⎢ ⎢ − 1 ⋱ − 1 ⎤ ⎦ ⎥ ⎥ ‚ B 1 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ D 1 W 1 ⋮ W i − 1 ⋮ W n − 1 W n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ‚ B 2 = ⎡ ⎣ ⎢ ⎢ G 1 ⋮ G n − 1 ⎤ ⎦ ⎥ ⎥ ‚ C 1 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ − F 0 c 0 β 1 η 1 α 2 η 2 E 1 − F 0 c 0 β 2 η 2 α 3 η 3 E 2 − F 0 c 0 β 3 η 3 ⋮ α i η i E i − 1 − F 0 c 0 β i η i ⋮ α n η n E n − 1 − F 0 c 0 β n η n W ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ , A 1 = [ α 1 η 1 - 1 0 α 2 η 2 - c * β 2 η 2 - 1 0 - c * β 3 η 3 α 3 η 3 - c * β 3 η 3 - 1 ⋮ ⋮ ⋮ ⋮ ⋱ 0 - c * β i η i ⋯ - c * β i η i α i η i - c * β i η i - 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 0 - c * β n η n ⋯ ⋯ - c * β n η n - c * β n η n α n η n c * β n η n - 1 0 1 ⋯ 1 ⋯ 1 1 1 ] , A 2 = [ c * ω 1 0 0 c * ω 2 c * ω 2 0 0 ⋮ ⋱ 0 0 c * ω i ⋯ ⋯ c * ω i 0 0 ⋮ 0 0 c * ω n - 1 ⋯ ⋯ ⋯ ⋯ c * ω n - 1 0 0 ] ‚ A 3 = [ 0 α 2 η 2 α 3 η 3 ⋱ α i η i ⋱ α n η n ] ‚ A 4 = [ - 1 ⋱ - 1 ] ‚ B 1 = [ D 1 W 1 ⋮ W i - 1 ⋮ W n - 1 W n ] ‚ B 2 = [ G 1 ⋮ G n - 1 ] ‚ C 1 = [ - F 0 c 0 β 1 η 1 α 2 η 2 E 1 - F 0 c 0 β 2 η 2 α 3 η 3 E 2 - F 0 c 0 β 3 η 3 ⋮ α i η i E i - 1 - F 0 c 0 β i η i ⋮ α n η n E n - 1 - F 0 c 0 β n η n W ] ,
C 2 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ 0 c ∗ ω 2 E 1 ⋮ c ∗ ω n − 1 ( E 1 + E 2 + ⋯ + E n − 2 ) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ C 2 = [ 0 c * ω 2 E 1 ⋮ c * ω n - 1 ( E 1 + E 2 + ⋯ + E n - 2 ) ]
各分块矩阵的阶数分别是:A 1 为n +1阶方阵, A 2 为 (n -1) × (n +1) 阶矩阵, A 3 为 (n +1) × (n -1) 阶矩阵, A 4 为 (n -1) 阶方阵;B 1 和C 1 均为 (n +1) ×1阶矩阵;B 2 和C 2 均为 (n -1) ×1阶矩阵。
《2.2 新算法简介》
2.2 新算法简介
解矩阵方程需给系数矩阵 (A 1 、A 2 和A 3 ) 中的元素和列向量 (C 1 和C 2 ) 中的元素赋初值, 这些元素与α i 、β i \, ω i \, η i 等量有关, 这些量又与t ′i 、t i 、H i 、w i 等未知量有关, 因而矩阵方程实质是非线性方程组。若给定α i 、β i 、ω i 、η i (η i 与w i 有关时给定w i ) 等变量的值, 则方程组转化为线性方程组。本文提出一种新算法——迭代法结合解线性方程组的矩阵法求解模型, 即先凭经验确定有关变量的初值, 将问题视为线性方程组;用高斯-约旦消去法[5 ] 解线性方程组得出D 1 、W i 、G i ;再按等面积原则得出Δt i [由式 (9) 可知计算Δt i 也必须给r i 和Δt 赋初值];最后用传热速率方程计算A i , 并判断是否符合精度要求, 若不符合则用修正后的有关变量值进行下一轮迭代计算直至收敛。新算法的计算框图如图2所示。
《图3》
图2 新算法程序框图
Fig.2 The new algorithm flow chart
据上述分析, 用迭代法结合矩阵法求解必须解决α i 、β i 、ω i 、w i 、r i 、Δt 的初值如何确定。根据文献[1 ] 的介绍, α i 近似为1, β i 的范围一般为0.0025~0.025 (第1效为沸点进料则β 1 =0) , 因而α i 、β i 的初值不难确定;ω i 的初值可参照β i 的范围取;r i 的初值可取2250×103 ;w i 的初值可按等蒸发量原则[1 ] 由式 (1) 和式 (2) 确定;Δt 的初值可取为t ′0 -t ′K -3。综上所述, 新算法所需的初值很容易确定。
《2.3 新算法的收敛性》
2.3 新算法的收敛性
迄今为止, 有关复杂并流多效蒸发系统的数学模型及求解方法未见报道, 为了检验新算法的优越性, 令块矩阵C 1 中的E i 为零, 并将块矩阵A 2 、A 3 、A 4 、B 2 和C 2 令为零矩阵, 则矩阵方程简化为没有冷凝水闪蒸和额外蒸汽引出的常规并流多效蒸发的矩阵。用新算法解文献[3 ] 报道的常规三效并流蒸发系统, 经过8次迭代即收敛, 而文献[3 ] 采用Broyden法求解需要赋12个变量的初值 (前已述及, 该法有些初值是难以确定的, 因而会影响迭代的稳定性和收敛性) , 经过19次迭代才收敛。由此可见, 新算法与目前常用的Broyden法相比具有对初值要求不高、收敛稳定性好、收敛速度快等优点。特别是当效数多、数学模型为大型方程组以及在多效蒸发系统的优化设计时, 新算法的优越性将更为明显。
《3 计算实例》
3 计算实例
兹有一并流加料的三效蒸发系统, 将含蔗糖w =10%的水溶液浓缩到50%, 沸点升高用下式Δ′i =1.78w i +6.22w 2 i i 2 (℃) 估算。加热用的饱和蒸汽压力为205 kPa (t ′0 为121.1 ℃) , 冷凝器的压力为13.7 kPa (t ′K 为51.3 ℃) 。进料速率为6.3 kg·s-1 , 温度为26.7 ℃, 溶液比热容为c 0 =4190-2350w 0 J·kg-1 ·K-1 。估算出的传热系数分别为:K 1 =3123 W·m-1 ·℃-1 , K 2 =1987 W·m-2 ·℃-1 , K 3 =1136 W·m-2 ·℃-1 。上述设计参数取自文献[4 ] , 文献[4 ] 的算例不考虑Δ″i 的影响并取η i =1, 本文按蒸发器内溶液液面高度为1.5 m计算Δ″i , 并取η i =0.98。假设各效的传热面积相等, 分别计算下列操作条件下的传热面积、生蒸汽的消耗量:1) 没有冷凝水闪蒸;2) 有冷凝水闪蒸;3) 假设用换热器将物料从26.7 ℃预热到90 ℃预热所需的热量为1577.2 kW, 若完全用生蒸汽加热要消耗加热蒸汽D 0 = 0.7159 kg·s-1 , 此部分热量由以下几个途径获得:①直接用生蒸汽作为加热源;②取第一效的额外蒸汽作为加热源;③取等量的一效和二效的额外蒸汽作为加热源。①、②、③三种情况蒸发系统均有冷凝水闪蒸。
用迭代法结合矩阵法求解, 将已知设计参数输入程序并运行得到模拟计算结果如表1所示。
从表1可知:对本例三效并流蒸发系统, 采用冷凝水闪蒸 (序号2) 比无冷凝水闪蒸 (序号1) 可节省加热蒸汽消耗量3.41%, 而传热面积仅增加1.44%。由于多效蒸发系统中加热蒸汽费用占年总费用的主要部分 (对于蒸发器材质为不锈钢或碳钢时约占88%~96%[6 ] , 故采用冷凝水闪蒸是节省年费用的有效措施;若同时采用冷凝水闪蒸和引出额外蒸汽预热原料液, 则节能效果更为显著, 以序号3中的③情况为例, 可节省加热蒸汽消耗量16.62%, 而各效传热面积还减小11.29%。
表1 算例模拟计算结果
Table 1 Calculation results of cocurrent three effect evaporation
《图4》
《4 结论》
4 结论
1) 建立了冷凝水闪蒸和额外蒸汽引出的复杂并流多效蒸发系统的数学模型, 常规并流多效蒸发系统只是该模型的一个特例。
2) 迭代法结合矩阵法是求解多效蒸发系统问题的一种新的有效方法, 与目前常用的Broyden法相比, 具有对初值要求不高、收敛稳定性好、收敛速度快等优点。
3) 本文的模型能够定量计算冷凝水闪蒸和额外蒸汽引出对多效蒸发加热蒸汽消耗量的影响。算例表明在三效并流蒸发时采用冷凝水闪蒸和额外蒸汽引出预热原料液可节省加热蒸汽消耗量16% 左右, 因而在并流多效蒸发系统中, 除末效以外的各效均应考虑冷凝水闪蒸和引出额外蒸汽预热原料液。