《1 引言》

1 引言

在化工过程的管路计算中, 通常需要知道摩擦因数并由Fanning公式计算流体流动过程中的机械能损失, 从而为流体输送机械的选型提供依据。工程上流体在粗糙管内的流动大多数是湍流流动, 其摩擦因数λ是雷诺数Re和相对粗糙度e/d的二元非线性函数。Colebrook方程 [1]:

《图1》

 

是得到工程界普遍认可、适用范围广 (Re=4×103~108, e/d=5×10-2~10-6) 、精度高 (与Murin等人实验数据的平均相对误差在5%以内 [2]) 的关于λRee/d之间的函数关系表达式, 但该方程是隐式方程, 计算摩擦因数要用迭代或试差的方法求解, 很不方便。Moody把Colebrook方程绘制成Moody摩擦因数图 [3], 该图直观明了, 使用方便, 但查图数据误差大。

随着计算机辅助设计CAD技术在化工领域的广泛应用, 用方程的形式来表达过程状态变量随状态条件的连续变化显得更为重要。为了克服Colbrook方程使用不方便的缺点, 许多研究者提出了多种形式的摩擦因数显式方程 [2,4,5,6,7,8,9,10,11,12,13,14]。在上述方程中, Altshul方程 [5]:

《图2》

 

形式最简单。对光滑管 (相对粗糙度e/d=0) , 该方程简化为著名的Blasius方程 [15]:

《图3》

 

但Altshul方程对Colebrook方程的总体平均误差高达8.2% [13]。文献[15]将Altshul方程中的指数0.25改为0.23、系数0.11改为0.1, 扩大了该方程的使用范围, 并得到了该方程适用于光滑管及阻力平方区 (λRe无关) 的两个简化式。其余方程 [2,4,6,7,8,9,10,11,12,13,14]均保留了Colebrook方程中的对数项, 形式复杂。李春喜和黄大铿方程 [13]:

《图4》

 

是其中形式相对简单、精度高 (与Colebrook方程的总平均相对偏差为0.07% [13]) 的一个方程。Romeo等人的方程 [14]:

《图5》

 

是其中形式最复杂、精度最高 (与Colebrook方程的平均相对偏差小于0.05% [14]) 。但这些方程均不易于简化成适用于光滑管或阻力平方区的方程。

笔者对先前提出的多元非线性多项式智能拟合法稍加改进 [16], 将智能拟合法应用于拟合Colebrook方程解的结果, 希望得到一个形式简单、精度高, 且便于简化成适用于光滑管或阻力平方区的计算管内湍动流体摩擦因数的显式新方程。

《2 智能拟合法》

2 智能拟合法

理论上, 任何复杂的多元非线性函数均可用拟合式来逼近。然而, 多元非线性拟合时如何筛选形式简单、精度高的最适宜拟合函数表达式是一项重要但又困难的工作。目前常用的多元非线性多项式拟合方法主要有 [17]:高斯-牛顿法、正交多项式回归法、阻尼最小二乘法、最速下降法、共轭梯度法、变尺度法、单纯型法等。这些常用的拟合方法存在以下两个问题:1) 它们大部分都是对固定形式的多项式的拟合, 多项式的项数及形式均由人工凭经验选定, 要找到最适宜的拟合函数表达式是及其困难的;2) 它们大多是对i=1nxi(i=012n)形常规幂多项式的拟合, 采用这种单一形式多项式拟合往往需要很多项才能达到要求的精度, 特别对于多元非线性多项式更是如此, 项数将会随着方次的升高成几何级数增加。如:二元多项式i=0n1j=0n2aijxiyj, 如果最高次n1, n2都取8, 拟合出来的项数就会有81个, 非常庞大, 不便使用。笔者在文献[16]中提出的一种多元非线性多项式新拟合法——智能拟合法, 可以较好地解决传统的多元非线性拟合方法存在的问题。

《2.1二元非线性多项式线性化处理及参数求解》

2.1二元非线性多项式线性化处理及参数求解

摩擦因数λ是雷诺数Re和相对粗糙度e/d的二元非线性函数, 在拟合前这种函数的适宜表达式 (包括项数、项形式、项系数和指数) 是未知的。为了拟合得到这种非线性函数的最适宜表达式, 首先必须构造非线性多项式的形式。为此, 令x=e/d, y=Re。考虑的常数项及具有整型次方的项的形式主要有:b0i=1n1j=1n2aijxiyji=1n3aixii=1n4aiyii=1n5j=1n6aij[lg(x)]i[lg(y)]ji=1n24aiy-i等, 共有16种形式的项 (其余未列出的项形式详见文献[16]) 。笔者另外考虑了具有实型次方的项如xd1, yd2, xd3yd4, (x+d5/y) d6等共4种形式。上述两大类型项共有20种项形式, 各项中x, y为待拟合函数的自变量, b0是常数项, nk (k=1, 2, 3, …, 24) 为各种形式项的整型最高次方, dl (l=1, 2, …, 6) 为实型参数。其次, 在给定项数、项形式、项指数时须求出项系数。矩阵法是拟合多元线性多项式的常用算法 [17], 它用最小二乘法推导得出具有最小误差的多项式拟合标准矩阵, 通过解矩阵直接求出误差最小的多项式各项系数, 没有其他一些方法的复杂递推或迭代过程, 算法简单高效。由于构造的项形式是二元非线性多项式, 为求出项系数, 将二元非线形多项式作线性化处理后亦可用矩阵法求解二元非线性多项式拟合问题 [16]。把前述的二元非线性项看成一个变量和一个常数的乘积, 即可转化成多元线型多项式。如对于二元高次多项式i=1n1j=1n2aijxiyj, 令b (i-1) ×n2+j=aijX (i-1) ×n2+j=xiyj, 即将其转化为线性项, 对其它形式的非线性项亦可作类似的处理将它们转化为线性项。再令Y=λ, 则摩擦因数λ的二元非线性方程可写成如下的n元线性方程:

《图6》

 

用最小二乘法可求得方程 (6) 的标准矩阵 [16,17]

《图7》

 

将解Colebrook方程得到的λ, Re, e/d之间的数据处理成m组数据XijYi (j=1, 2, …, n, 代表自变量代号;i=0, 1, 2, …, m, 代表数据组号) 代入标准矩阵方程 (7) 并解之, 即可得到要拟合的项系数bj (j=0, 1, 2, …, n) 。

上述方法只能求出二元非线性多项式在给定项数、项形式、项指数时与待拟合数据之间“误差最小”的项系数, 此“误差最小”只是局部最优, 并非全局最优。要得到全局最优, 即得到与待拟合数据之间误差最小的最适宜拟合函数表达式 (包括项数、项形式、项指数和系数) , 必须用下面介绍的遗传算法实现。

《2.2用遗传算法实现智能拟合》

2.2用遗传算法实现智能拟合

由上述模型可知, 决定拟合函数表达式的项数、项形式和项指数的是24个整型变量 (n1, n2, …, n24) 和6个实型变量 (d1, d2, …, d6) , 共30个变量, 称之为控制变量。若能用一种最优化方法搜索这30个控制变量的值并结合矩阵法求项系数, 根据拟合式与待拟合数据之间的误差大小决定项的取舍, 最终就能搜索得到全局最优的最适宜拟合函数表达式, 实现智能拟合。然而, 面对如此庞大的变量阵容, 用普通的优化算法是很难实现智能拟合的。遗传算法是一种新兴的仿生算法, 具有快速寻优和全局优化的特点, 笔者在文献[16]中用杰出者选择遗传算法实现了具有整数方次的二元非线性多项式的智能拟合 (该拟合问题的控制变量是前述的24个整型变量) , 效果很好, 故作者在整数次方项的搜索和取舍问题以及选择、交叉、变异、整型二进制编码等策略方面均沿用文献[16]中的方法, 下面着重说明具有实数方次项的搜索、取舍及实型变量的编码等问题。

对具有实数次方的项, 一方面要搜索项内参数, 另一方面要控制项的取舍问题, 不能像文献[16]中控制整数次方项那样处理, 对于每一个实参项类型都要有一个对应的整型变量来专门控制项的取舍。项的取舍须考虑取舍概率问题。舍弃率太大会导致算法不稳定, 收敛困难;舍弃率太小会使得该项成为整个模型的累赘, 容易陷入局部最优。整数次方的项的取舍变量是整型参数 (最高方次nk) 本身, 用3位二进制表示, 值取0时舍弃该项, 舍弃率为1/8。为与整型项协调统一和方便控制, 具有实参的项的舍弃率也取1/8较好, 实参项的控制变量取[0, 7]范围的整数 (用3位二进制编码) , 值取0时舍弃该项。4个具有实型次方的项形式, 得加4个整型控制变量 (v1, v2, v3, v4) 来控制实参项的取舍问题。

实型变量和整型变量的编码有所不同, 也采用二进制编码的方法, 每个变量用20位二进制表示, 用6×20的二维数组表示这6个实型变量的代码。变量与代码的对应关系可用式 (8) 表示。

《图8》

 

式中, di为实型变量;Si, j是二进制代码, 0或1;i为变量序号, ki=1, 2, …, 6;j是代码位号, j=1, 2, …, 20。aibi为变量的取值范围, 其中ai<di<bi

初始种群一般随机产生, 但产生的概率可根据变量特征人为选定。根据变量特征, 取1和0的产生概率为1∶3。根据经验, 种群规模取80个, 其中每代新产生79个个体, 1个是上一代保留下来的最佳个体。文献[16]中采用均匀交叉法, 均匀交叉法易于保持个体的多样性, 对于大量短代码变量的寻优, 处理简单、效果好。单点交叉每对变量的交叉点只有一个, 易于保留父个体的优秀模式, 有利于长代码的交叉。采用2种交叉法相结合的方法, 整型变量对应的代码采用随机交叉, 交叉率取1/3;实型变量对应的代码采用单点交叉法, 交叉概率取3/4。变异概率是控制进化历程的重要参数, 太小会使算法过早收敛而出现早熟现象, 太大则降低算法效率, 收敛慢, 根据经验取0.02。

《3 管内湍动流体摩擦因数的显式方程》

3 管内湍动流体摩擦因数的显式方程

从Colebrook方程适用范围内 (Re=4×103~108, e/d=5×10-2~10-6) 均匀选取400个点用Colebrook方程算得摩擦因数值, 应用上述智能拟合法进行拟合, 得到摩擦因数计算的显式新方程:

《图9》

 

适用于Moody摩擦系数图 (该图是根据Colebrook方程绘制的 [3]) 4×103Re≤1×108的整个范围, 包括湍流区、完全湍流区、光滑管 (e/d=0) ;适用的相对粗糙度范围为:0≤e/d≤0.05。式 (9) 对Moody图的上述范围的平均误差为0.5%, 最大误差小于2%。 对于光滑管 (e/d=0) , 式 (9) 可简化为式 (10) ;对于完全湍流区 (Reλ影响可以忽略) , 式 (9) 可简化为式 (11) 。

《图10》

 

将文献[2]摩擦系数λ的实验值与各公式计算的λ值对照得表1, 对其进行误差分析见表2。由表2可知, λ的计算值对实验值的偏差, 式 (2) 最大, 式 (1) , 式 (4) 和式 (5) 居中, 新方程式 (9) 最小。由表2还可知, 对Colebrook方程的偏差, 式 (2) 最大, 但其形式最简单;式 (5) 最小, 但其形式最复杂;式 (4) 也很小, 其形式亦很复杂;式 (9) 较小, 其形式亦较简单。由于Colebrook方程是个半经验式, 对实验数据存在较大偏差, 所以式 (4) , 式 (5) 对其偏差虽然极小意义却不大, 且式 (4) , 式 (5) 形式复杂, 使用不便, 不能简化为完全湍流区及光滑管的计算公式。式 (2) 虽然形式简单, 但其对实验数据及Colebrook方程的偏差很大, 精度低, 不能满足工程计算的要求。式 (9) 对实验数据及Colebrook方程的偏差均很小, 精度高, 能满足工程计算的要求, 且其形式简单、使用方便、易于简化成光滑管的式 (10) 及完全湍流区的式 (11) , 故式 (9) 是一个精度高、形式简单的计算管内湍动流体摩擦因数的显式新方程。


  

表1 有关公式λ的计算值与实验值对照表  

Table 1 Comparison between the experimental data and the results calculated by formulas for λ

 

《图11》

表1 有关公式λ的计算值与实验值对照表


  

表2 有关λ公式计算值对实验值及Colebrook方程计算值的误差分析表  

Table 2 Deviation of the formulas for λ to the experimental data and the Colebrook equation

 

《图12》

表2 有关λ公式计算值对实验值及Colebrook方程计算值的误差分析表

《4 结论》

4 结论

1) 智能拟合法是一种新的非线性拟合法, 具有可以自动搜索最适宜拟合函数表达式、拟合精度高、拟合速度快等优点, 可以广泛应用于解决科学与工程领域的众多二元非线性拟合问题;

2) 与前人获得计算摩擦因数显式方程的途径不同, 笔者应用智能拟合法拟合Colebrook隐式方程得到形式简单、精度高、适用范围广的计算管内湍动流体摩擦因数的显式新方程, 使用方便;

3) 新方程对实验数据的相对偏差为2.3%, 可与Colebrook方程相媲美, 完全满足工程设计计算需要, 对Colebrook方程的偏差为0.5%, 完全可以作为Colebrook方程代替公式;