公众对环保问题的重视与新环保法的实施，使燃油标准的提升成为必然趋势，从而对炼油过程提出了开展分子水平精细化管理的要求。从分子水平认识石油馏分是建立分子尺度油品转化过程模型的前提，也是实现炼油过程分子管理的基础。当前针对石油馏分的分子层次表征方法主要有以下3种：
1) 同分异构分子重构模型(MTHS)：模型由Peng[1] 在1999年提出，将石油馏分按照碳数、结构特征(同分异构体)进行分类，使用矩阵来表述各组分，矩阵的行和列分别代表碳数和特征结构，矩阵的每个元素视为一个虚拟组分。Gomez-Prado等[2]利用MTHS对重油馏分进行了分子重构，模拟计算结果与实测值吻合较好。但该模型需要预先建立虚拟组分库，如果组分选择不当，将会对计算的精度产生很大影响。
2) 统计重构模型：该模型利用特定属性或官能团(芳环数、侧链数、侧链长度、环烷环数等)区分石油馏分中的分子，并利用概率密度函数表达组分分布。Hudebine等[3]对该模型进行了完善和扩展，提出了最大熵值统计重构模型(REM)，利用虚拟组分的信息熵最大化原则，使模拟混合物的微观特性也趋近所要表征的石油馏分。彭辉等[4]利用此方法对乙烯裂解原料进行分子重构，产生的等效分子的常规物性计算值与标定值吻合较好。该模型考虑了宏观性质和微观特性，也达到了不错的计算精度，但不易表征含杂原子较多的石油馏分。
3) 结构导向集总模型(SOL)：该模型用向量来表示分子，向量的元素代表了分子中结构的个数。这种方法是分子尺度上的集总，这样利于降低模型的复杂性，将体系中的分子由数以百万计减少到几千。虚拟分子不但包括烷烃、环烷烃和芳烃结构，还包括作为中间产物或一次反应产物的烯烃及环烯烃，同时考虑了含硫、氮等杂原子的化合物。Ghosh和Jaffe[5, 6]利用SOL模型估算了203种柴油的十六烷值和1 471种汽油的辛烷值，结果表明预测误差均小于2%。该模型利用随机生成分子的形式，不需要预先建立分子库，同时基本结构中有若干杂原子结构，非常适用于含杂原子较多的石油馏分的表征。
本研究采用基于结构导向集总模型的方法，选取22个特征结构表征石油馏分，通过Monte Carlo模拟产生虚拟分子，同时借鉴统计重构模型中的信息熵值最大化思想，用优化的方法调整虚拟分子的含量，使得重构分子不仅在宏观性质上反映实际原料特性，在微观特性上也趋近所表征原料。
1 石油馏分宏观性质石油馏分种类很多，此处选取减压蜡油(VGO)进行表征。表 1为2种减压蜡油的宏观性质，其中各元素、残炭和族组成的百分数均指质量分数。
| 石油馏分 | 密度ρ20 | 运动粘度ν80 | 平均分子量M | 残炭质量分数 | 元素组成及其质量分数/% | 族组成及其质量分数/% | |||||
| g·cm-3 | mm2·s-1 | % | C | H | S | N | 链烃 | 环烷烃 | 芳烃 | ||
| 大连VGO | 0.922 6 | 22.98 | 407.1 | 0.30 | 86.91 | 12.33 | 0.54 | 0.17 | 9.9 | 47.1 | 43.0 |
| 苏丹VGO | 0.926 1 | 15.21 | 346.0 | 0.04 | 86.93 | 12.37 | 0.26 | 0.19 | 12.7 | 48.6 | 38.7 |
结构导向集总方法构建虚拟分子的核心是选用一些结构来表达分子中的关键模块，然后使用这些结构拼装分子。构建的分子用向量表示，向量的分量是各个特征结构的数量。Jaffe[7]认为，这样模拟得到的虚拟分子有可能在真实的原料油中并不存在，但是大量虚拟分子的宏观性质表现出与真实原料的高度相似性，证明这种方法是正确可行的。
2.1 虚拟分子结构的选择VGO中除了碳、氢这2种主要元素外，还含硫、氮、氧等杂原子和重金属。根据表 1中的标定数据，可以对问题进行适当的简化，只考虑碳、氢、硫、氮和氧这5种元素。最终选取包含以上5种元素的22个结构表达虚拟分子[7]，各结构的表示符号以及结构中各个原子的计量个数表 2所示。
| 原子种类 | A6 | A4 | A2 | N6 | N5 | N4 | N3 | N2 | N1 | R | br | me | IH | AA | NS | RS | AN | NN | RN | NO | RO | KO |
| C | 6 | 4 | 2 | 6 | 5 | 4 | 3 | 2 | 1 | 1 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | -1 | 0 | -1 | 0 | 0 |
| H | 6 | 2 | 0 | 12 | 10 | 6 | 4 | 2 | 1 | 2 | 0 | 0 | 2 | -2 | -2 | 0 | -1 | -1 | 1 | -2 | 0 | 0 |
| S | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| N | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 |
| O | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 |
A6：表示芳环的主要成环结构，可以单独存在，当其他分量为0时代表一个单体苯。
A4：表示含4个碳的芳环，不能单独存在，只能连接在A6或另一个A4结构上。
A2：表示含2个碳的芳环，不能单独存在，只能依附于A6和A4，用来连接相邻的缩合芳环，是形成稠环芳烃的主要部分。
2.1.2 脂环结构N6和N5：脂环的主要成环结构，可以单独存在，当其他分量为0时分别代表环己烷和环戊烷。
N4、N3、N2、N1：分别表示含4、3、2、1个碳的脂环，不能单独存在，在分子结构中连接在芳环或者脂环上的附属脂环。
2.1.3 烷基结构R：表示除芳环和脂环之外的烷基总碳数。
br：表示烷基侧链的个数。
me：表示直接连在环结构上的甲基数目，但R=1时，me=0。
2.1.4 辅助结构IH：表示不饱和度(芳环除外)，每增加一单位IH，分子中减少2个氢。无脂环时，IH为1代表链烷烃，IH为0代表单烯，IH为-1代表双烯；有脂环时，IH为0代表脂环，IH为-1代表环烯。
AA：A6、N6和N5之间的桥键。
2.1.5 杂原子结构NS、NN、NO：表示脂环或烷烃上的硫、氮和氧原子。
RS、RN、RO：表示碳氢原子间的硫、氮和氧原子。
AN：表示芳环上的氮原子。
KO：表示羰基或醛基上的氧原子。
2.2 各结构存在条件、约束条件以及取值范围的规定从上述结构描述中可以看出，并不是所有结构向量都能表示有意义的分子，各结构的数量需要符合一定的条件，结合相关资料[8, 9]中对VGO性质的描述，本文作如下规定：
A6：最小值为0，最大值为3。
A4：当A6≥1时存在，最小值为0，最大值为3。
A2：当A4≥2时存在，最小值为0，最大值为A4-1。
同时，还规定A6+A4+A2≤5，即分子中芳环总数不大于5。
N6、N5：最小值均为0，最大值均为3。
N4：当A6+N6≥1时存在，最小值为0，最大值为3。
N3：当A6+N6+N5≥1时存在，最小值为0，最大值为3。
N2：当A4+N4+N3≥1时存在，最小值为0，最大值为A4+N4+N3。
N1：当A4+N4+N3≥2时存在，最小值为0，最大值为A4+N4+N3-1。
同时，还规定N6+N5≤3，即分子中单独的五六元环总数不大于3；N6+N5+N4+N3+N2+N1≤8，即分子中脂环总数不大于8。
R：当A6+N6+N5=0(无环结构)时，最小值为6，最大值为40；当A6+N6+N5>0(有环结构)且环结构碳数≥40时，R=0；当A6+N6+N5>0(有环结构)且环结构碳数 < 40时，最小值为6，最大值为40-环结构碳数。
br：当R≤3时，br=0；当R=4时，最小值为0，最大值为1；当R>4时，最小值为0，最大值为R-4和4中的小值，即min(R-4，4)。
me：当A6+N6+N5>0时存在。当R≤1时，me=0；当R=2时，最小值为0，最大值为1；当R=3 时，最小值为0，最大值为2；当R≥4时，若br=0则最小值为0，最大值为min(R-1，4)，若br≥1则最小值为0，最大值为min(R-br-3，4)。
IH：当A6+N6+N5>0时为0；当A6+N6+N5=0时为1。
AA：当A6+N6+N5≥2(有2个以上独立的环结构)时存在，AA= A6+N6+N5-1。
NS、NN：当A6≥2时存在，最小值均为0，最大值为均为min(A6-1，AA)。
RS、RN：最小值为0，最大值为3。
AN：当A6≥1时存在，最小值为0，最大值为A6+A4+A2。
NO、RO、KO：最小值均为0，最大值均为1。
同时，还规定单个分子中硫、氮、氧原子不同时出现，且除羧酸外，杂原子总数不大于1。
各结构的相关规定是以实际原料为依据进行的，可以根据实际情况进行相应调整。
2.3 虚拟分子Monte Carlo抽样一般认为在22个基本结构中，除IH和AA这2个辅助结构外，其余20个结构在每个分子中的数量均服从χ2分布且相互独立。χ2分布的概率密度函数为：
| $p\left( x \right) = \frac{{{{\left( {x - {x_{\min }}} \right)}^{\frac{r}{2} - 1}} \cdot {{\rm{e}}^{ - \frac{{x - {x_{\min }}}}{2}}}}}{{\Gamma \left( {\frac{r}{2}} \right) \cdot {2^{\frac{r}{2}}}}};\;\;\;r = \mu = \frac{{{\sigma ^2}}}{2}.$ | (1) |
其中Γ(α)为Gamma函数，即
| $\Gamma \left( \alpha \right) = \int_0^\infty {{t^{\alpha - 1}}{{\rm{e}}^{ - t}}{\rm{d}}{\rm{.}}} $ | (2) |
令式(3)成立，有
| $\int_{{x_{\min }}}^{{x_{\max }}} {\frac{{{{\left( {x - {x_{\min }}} \right)}^{\frac{r}{2} - 1}} \cdot {{\rm{e}}^{ - \frac{{x - {x_{\min }}}}{2}}}}}{{\Gamma \left( {\frac{r}{2}} \right) \cdot {2^{\frac{r}{2}}}}}} {\rm{d}}x = Q.$ | (3) |
则可以对每个分子中各个结构的数量进行抽样，抽样公式如式(4)所示，即
| $\int_{{x_{\min }}}^x {\frac{{{{\left( {x - {x_{\min }}} \right)}^{\frac{r}{2} - 1}} \cdot {{\rm{e}}^{ - \frac{{x - {x_{\min }}}}{2}}}}}{{\Gamma \left( {\frac{r}{2}} \right) \cdot {2^{\frac{r}{2}}}}}} {\rm{d}}x = {\rm{RN}}.$ | (4) |
其中：x为每次抽样得到的样本值；xmin、xmax分别为结构数量的最小值、最大值；r为x的数学期望；μ为χ2分布的自由度；RN是随机数，服从[0,Q]内的均匀分布。
为了确保得到的抽样值在设定的最大值和最小值之间，每个分子中各结构数量的具体抽取方法如下：首先给定各结构所服从的χ2分布的参数(自由度)；然后根据式(3)计算出各结构数量最大值和最小值的分布函数值之差Q；最后随机抽取RN，根据式(4)进行逆χ2分布计算得到抽样值x。上述过程每循环一次即可得到一个结构的数量，用22个结构的数量便可表示一个虚拟分子，如此循环抽样N次则生成N个虚拟分子。
由于实际生产中往往难以得到VGO的所有性质，因此对比性质应从可获得的标定数据中选择，如选取以下8项：C、H、S、N的质量分数，芳烃、环烃、链烃的质量分数以及平均分子量。
抽样得到若干虚拟分子之后，由于已知每个分子的结构特性，因此可以计算出每个分子的上述性质，从而得到全部虚拟分子的混合性质。相关性质计算值与标定值的相对偏差平方和为：
| $\begin{array}{*{20}{c}} {y = {{\sum\limits_{i = 1}^4 {\left( {\frac{{{\rm{EL}}{{\rm{E}}_{i,{\rm{cal}}}} - {\rm{EL}}{{\rm{E}}_{i,\exp }}}}{{{\rm{EL}}{{\rm{E}}_{i,\exp }}}}} \right)} }^2} + }\\ {{{\sum\limits_{i = 1}^3 {\left( {\frac{{{\rm{PION}}{{\rm{A}}_{i,{\rm{cal}}}} - {\rm{PION}}{{\rm{A}}_{i,\exp }}}}{{{\rm{PION}}{{\rm{A}}_{i,\exp }}}}} \right)} }^2} + }\\ {{{\left( {\frac{{{M_{{\rm{cal}}}} - {M_{\exp }}}}{{{M_{\exp }}}}} \right)}^2}.} \end{array}$ | (5) |
其中：脚标cal和exp分别表示计算值和标定值；ELE表示元素(C、H、S、N)的质量分数；PIONA表示各族组成(链烃、环烃、芳烃)的质量分数；M表示平均分子量。因为是初步估算分子结构以及组成，可通过优化r使y小于阈值δ，从而得到大致能表现实际原料宏观特性的虚拟分子。优化方法采用模拟退火法，基于Matlab编程实现，最终得到各结构对应的r值见表 3和表 4。
| 结构 | r | 结构 | r |
| A6 | 0.850 3 | br | 1.324 2 |
| A4 | 0.858 8 | me | 1.329 6 |
| A2 | 0.850 6 | NS | 0.182 6 |
| N6 | 1.211 7 | RS | 0.185 8 |
| N5 | 1.207 5 | AN | 0.110 0 |
| N4 | 1.211 8 | NN | 0.105 2 |
| N3 | 1.213 6 | RN | 0.104 1 |
| N2 | 1.216 7 | NO | 0.131 9 |
| N1 | 1.208 3 | RO | 0.131 0 |
| R | 16.701 0 | KO | 0.129 6 |
| 结构 | r | 结构 | r |
| A6 | 0.719 6 | br | 1.209 8 |
| A4 | 0.720 1 | me | 1.329 6 |
| A2 | 0.720 3 | NS | 0.083 5 |
| N6 | 1.115 8 | RS | 0.084 5 |
| N5 | 1.116 2 | AN | 0.113 6 |
| N4 | 1.120 0 | NN | 0.114 1 |
| N3 | 1.115 0 | RN | 0.114 0 |
| N2 | 1.115 7 | NO | 0.120 2 |
| N1 | 1.116 1 | RO | 0.118 9 |
| R | 13.200 0 | KO | 0.122 2 |
1948年，Claude Shannon[10]把熵的概念引入到信息论里，称为信息熵。信息熵是对不确定性的测量，通常用来表示系统的不确定度，在信息熵最大的时候随机变量的不确定性也越大，即更符合随机的情况。信息熵可用式(6)来表达，即
| $S\left( z \right) = - \sum\limits_{i = 1}^N {{z_i}\ln {z_i}} .$ | (6) |
其中：S表示信息熵；zi≥0,(i=1,2，…，N),$\sum\limits_{i = 1}^N { = 1} $.
在2.3节所提出的构造虚拟分子体系来表征原料油的模型中，虚拟分子是按照等摩尔处理的。根据信息熵最大化原理，最合理的分子组成，应使得信息熵最大。基于此，可在上节结果的基础上进一步优化，调整分子量。将预测的石油烃的分子组成zi的信息熵作为目标函数，宏观性质计算值与标定值偏差平方和小于阈值σj作为约束条件，即
| $\begin{array}{*{20}{c}} {\max S\left( z \right) = - \sum\limits_{i = 1}^N {{z_i}\ln {z_i}} ;}\\ {{\rm{s}}{\rm{.t}}{\rm{.}}{{\left( {{f_i} - \sum\limits_{i = 1}^N {{f_{i,j}}{z_i}} } \right)}^2} \le {\sigma _j},j = 1, \cdots ,J.} \end{array}$ | (7) |
其中：S(z)表示熵值；zi表示分子摩尔分数，0≤zi≤1(i=1,…,N)，$\sum\limits_{i = 1}^N {{z_i} = 1} $；fj表示可以获得的各常规物性数据，如元素质量分数、族组成质量分数等，j表示对比的性质数量。$\sum\limits_{i = 1}^N {{f_{i,j}}{z_i}} $表示通过重构产生的分子组成的总物性，而σj为预测物性的允许误差。相较于阈值δ，阈值σj要求各项性质计算值与标定值偏差更小。约束条件一般根据可获得的原料物性数据来决定。本文中使用的约束项如式(8)—式(10)所示，即
| ${\left( {E - \sum\limits_{i = 1}^N {{E_i}{z_i}} } \right)^2} \le {\sigma _1},\;\;\;\;E \in \left( {C,H,S,N} \right).$ | (8) |
| ${\left( {L - \frac{{\sum\limits_{i \in L} {{M_i}{z_i}} }}{{\sum\limits_{i = 1}^n {{M_i}{z_i}} }}} \right)^2} \le {\sigma _2},\;\;\;L \in \left( {P,N,A} \right).$ | (9) |
| ${\left( {M - \sum\limits_{i = 1}^N {{M_i}{z_i}} } \right)^2} \le {\sigma _3}.$ | (10) |
其中：E、L和M分别表示实际原料油的元素组成质量分数、结构族组成质量分数和平均分子量；Ei和Mi分别表示虚拟分子的元素质量分数和分子量。
这个优化问题是有约束非线性规划问题，可使用罚函数法将该有约束问题转化为无约束问题来处理，将约束条件加到目标函数中，新的目标函数如式(11)所示，即
| $\begin{array}{*{20}{c}} {F\left( {z,c} \right) = - S\left( z \right) + c \cdot g\left( z \right),}\\ {g\left( z \right) = \sum\limits_{j = 1}^J {{{\left( {{f_i} - \sum\limits_{i = 1}^N {{f_{i,j}}{z_i}} - {\sigma _j}} \right)}^2}} .} \end{array}$ | (11) |
其中c为惩罚因子。此步优化仍然选取模拟退火法，利用Matlab编程实现。基于信息熵理论的分子组成求解方法的框图如图 1，调整分子量后得到的原料的标定值与计算值的对比见表 5和表 6。从表 5和表 6可以看出得到了很好的计算结果。
|
| 图 1 基于信息熵理论的分子组成求解框图 |
| 结果 | 平均分子量M | 元素组成及其质量分数/% | 族组成及其质量分数/% | |||||
| C | H | S | N | 链烃 | 环烷烃 | 芳烃 | ||
| 标定值 | 407.10 | 86.91 | 12.33 | 0.54 | 0.17 | 9.90 | 47.10 | 43.00 |
| 计算值 | 403.03 | 86.68 | 12.30 | 0.54 | 0.17 | 9.97 | 47.91 | 42.12 |
| 相对误差/% | -1.01 | -0.27 | -0.24 | 0 | 0 | 0.70 | 1.69 | -2.09 |
| 结果 | 平均分子量M | 元素组成及其质量分数/% | 族组成及其质量分数/% | |||||
| C | H | S | N | 链烃 | 环烷烃 | 芳烃 | ||
| 标定值 | 346 | 86.93 | 12.37 | 0.26 | 0.19 | 12.7 | 48.6 | 38.7 |
| 计算值 | 348.74 | 86.79 | 12.43 | 0.26 | 0.19 | 12.35 | 49.01 | 38.64 |
| 相对误差/% | 0.79 | -0.16 | 0.49 | 0 | 0 | -2.76 | 0.82 | -0.16 |
至于抽取的分子数，由于样本量过小会影响模拟结果的准确性，而样本量过大会使程序运行的时间大大增加，因此需要进行综合考虑。现有研究[11-13]模拟的分子数基本在1 000~5 000个。
本研究分别在相同精度下模拟了1 000、2 000和 5 000 个分子，结果显示，随着分子数的增加，程序运行时间成倍增加，从几小时上升至几天。当虚拟分子数为 1 000 时，模拟得到的结果精度与程序运行时间可以得到较好的平衡。本文所列出的结果均为虚拟分子数为1 000时的运算结果。
4 结 论本文借鉴结构导向集总模型中对分子的表示方法，选取了包含杂原子硫、氮和氧的22个基本结构，使得分子的元素构成更加符合实际。同时使用Monte Carlo方法随机抽取虚拟分子，并利用信息熵最大化原则对各虚拟分子的质量分数进行了调整，使得模拟得到的分子不但在宏观性质上与原料很好地吻合，在微观特性上也趋于实际原料。计算值与标定值的相对误差不超过3%，说明此方法对于复杂反应体系的原料重构是可行的。
另外，本文通过优化分布参数和调整分子量，对2套不同的石油馏分进行了分子重构。计算得到的平均性质均能和标定数据很好地拟合，说明此方法能够广泛地适用于重构各种不同性质的复杂馏分原料。
| [1] | PENG Bin. Molecular Modeling of Petroleum Processes[D]. Manchester, UK:Institute of Science and Technology of the University of Manchester, 1999. |
| [2] | Gomez-Prado J, ZHANG Nan, Theodoropoulos C. Characterization of heavy petroleum fractions using modified molecular-type homologous series (MTHS) representation[J]. Energy, 2008, 33(6):974-987. |
| [3] | Hudebine D, Verstraete J J. Molecular reconstruction of LCO gasoils from overall petroleum analyses[J]. Chemical Engineering Science, 2004, 59P(22-23):4755-4763. |
| [4] | 彭辉, 张磊, 邱彤, 等. 乙烯裂解原料等效分子组成的预测方法[J]. 化工学报, 2011(12):3447-3451.PENG Hui, ZHANG Lei, QIU Tong, et al. Method of predicting equimolecular mixture of ethylene cracking feedstock[J]. CIESC Journal, 2011(12):3447-3451. (in Chinese) |
| [5] | Ghosh P, Hickey K J, Jaffe S B. Development of a detailed gasoline composition-based octane model[J]. Industrial & Engineering Chemistry Research, 2006, 45(1):337-345. |
| [6] | Ghosh P. Predicting the effect of cetane improvers on diesel fuels[J]. Energy & Fuels, 2008, 22(2):1073-1079. |
| [7] | Quann R J, Jaffe S B. Structure-oriented lumping:describing the chemistry of complex hydrocarbon mixtures[J]. Industrial & Engineering Chemistry Research, 1992, 31(11):2483-2497. |
| [8] | 韩崇仁. 加氢裂化工艺与工程[M]. 第一版. 北京:中国石化出版社, 2006.HAN Chongren. The Craft and Engineering of Hydrocracking[M]. 1st Ed. Beijing:China Petrochemical Press, 2006. |
| [9] | 方向晨. 加氢裂化[M]. 第一版. 北京:中国石化出版社, 2012. FANG Xiangchen. Hydrocracking[M]. 1st Ed. Beijing:China Petrochemical Press, 2012. |
| [10] | Shannon C E. A mathematical theory of communication[J]. Bell System Technical Journal, 1948, 27(3):379-423. |