离心泵蜗壳拓扑优化设计方法研究

徐胜利，张辉，金婷婷，杨成立，周演

(1. 大连理工大学能源与动力学院，辽宁 大连 116024； 2. 上海电气燃气轮机有限公司，上海 200240； 3. 大连理工大学工程力学系，辽宁 大连 116024)

蜗壳是离心泵的重要过流部件之一，其水力性能对离心泵的整体水力性影响显著.在离心泵输送流体过程中，蜗壳中的水力损失达到离心泵整体水力损失的一半以上[1]，蜗壳水力性能优化设计对提高离心泵的水力效率具有重要意义.传统的蜗壳一维设计方法简单方便，目前发展已经十分成熟，主要包括平均速度法、等环量法、阿基米德螺旋线方程法、结构方形法、不等距方形法等[2].这些一维设计方法需要基于一定的假设条件，如恒定角动量假设、恒定平均速度假设等[3].蜗壳由于受非轴对称几何形状的影响，其内部的瞬态流场呈非轴对称分布[4]，特别是受蜗壳隔舌的影响，使得蜗壳入口处的流动参数的实际分布与理想假设有差别[5].同时，在一些特殊设计条件下，如外部空间复杂、安装尺寸受限时，传统蜗壳设计方法难以满足对空间位置和尺寸限制的要求.因此，发展一种不依赖于已有设计经验，能够在外部空间受限条件下进行自动蜗壳优化设计的方法显得十分重要.

近年来，流体拓扑优化方法发展迅速.拓扑优化密度法通过多孔介质模型建立材料密度与流动阻力的关系，实现设计域中固体和流体之间的转化[6].起初，流体拓扑优化主要应用在低雷诺数流动问题的流道设计中[7-8].随着冻结湍流方法的提出[9]，湍流黏度的变化可以忽略，流体拓扑优化也能够应用到湍流问题的优化设计中[10-11].由于拓扑优化设计变量多，通常采用梯度优化算法求解优化问题，伴随方法是性能函数灵敏度计算的主要方法.伴随方法通过原始控制方程和性能函数导出伴随方程和伴随边界条件，通过求解伴随方程即可快速得到灵敏度的值.应用连续伴随法，已经实现了二维、三维管道流动的拓扑优化[10].近年来，流体拓扑优化方法被应用到离心泵叶轮的流道拓扑优化设计中[12-13].离心泵试验测试也验证了拓扑优化方法能够提升泵的整体性能[14].

文中采用流体拓扑优化方法进行离心泵的蜗壳流道优化设计，建立适用于外部空间受限条件下蜗壳设计的拓扑优化方法，基于开源软件OpenFOAM建立相应的优化问题分析和求解程序.针对3种不同的蜗壳出口位置，进行蜗壳流道优化设计，获得合理的蜗壳形状.

1 优化问题

1.1 优化列式

以离心泵的蜗壳为研究对象，构建了合适的设计域，并假设设计域为多孔介质模型，对其进行优化设计，优化问题可描述为

find:ρ=(ρ1,ρ2,ρ3,…,ρn),

s.t.:R[v,p,α(ρ)]=0,

(1)

0≤ρi≤1, (1≤i≤n),

式中：ρ为单元密度；下标n为设计域内单元总数目；J为目标函数，即设计域进出口的总压差与流体密度的比值，m5/s3；p为静压与流体密度的比值，m2/s2；v为流体速度；S为设计域的边界；n为边界面的单位法向量；R为流动控制方程；α为多孔介质孔隙度的倒数；V为单元体积；m为流体域所占设计域总体积的比例上限.

当单元密度ρ=1时，α=αmax，对应为固体单元；ρ=0时，α=αmin，对应为流体单元.对于ρ在0～1之间的中间密度值，通过基于Heaviside函数的非线性密度过滤方法，使中间密度值趋于0或1.流体在多孔介质中流动时受到的阻力被称为Darcy阻力，其大小与孔隙度的倒数α以及流体流动速度有关，可表示为f=αv.将Darcy阻力项加入不可压缩、稳态的Navier-Stokes方程中，即可得到文中所使用的控制方程

(2)

(3)

式中：ν为运动黏度，对于湍流，ν可认为是流体的运动黏度与湍流黏度之和.

1.2 密度过滤

在拓扑优化中，优化结果可能会出现棋盘分布现象，即固体单元ρ=1和流体单元ρ=0交替排列，类似于棋盘状图案，不能得到可用的优化设计方案.通过非线性密度过滤来控制设计域内的密度分布，能够消除优化过程中可能出现的棋盘分布现象.首先对单元某一特定范围内的全部单元密度进行加权平均，再以密度平均值作为范围中心这一单元的线性过滤密度.线性密度过滤如图1所示，图中C为线性密度过滤范围圆心，Rg为过滤半径.

图1 线性密度过滤示意图

(4)

式中：Xi为过滤范围内单元i中心点的坐标；XC为单元C中心点的坐标；w(Xi)为过滤范围内所有单元密度的权重，其与Xi的对应函数关系为

(5)

(6)

由于式(6)为分段函数，故需要对其进行连续化处理，如式(7)所示，以此保证灵敏度的连续性.

(7)

式中：β为控制函数凹凸程度的参数.图2为β取不同值的非线性函数曲线.

图2 连续化Heaviside函数

(8)

式中：αmin为多孔介质孔隙度倒数α的下限(文中为0)；αmax为α的上限，理论上αmax应趋近于无穷大，在实际计算中选取一个较大的数来近似.

1.3 伴随方程

采用梯度类优化算法求解优化问题，需要计算目标函数对所有设计变量的灵敏度

(9)

(10)

式中：Ω为整个设计域.拉格朗日乘子可表示为

(u,q)=(u1,u2,u3,q).

(11)

对拉格朗日函数关于设计变量求全导，可得

(12)

显然

(13)

目标函数J与设计变量ρ无关，因此式(13)中等号右端第1项为0.

目标函数J可以表示为设计域内部积分与边界积分之和，即

(14)

将方程(2)，(3)代入到式(13)中，结合分部积分公式，整理后可得

(15)

(16)

(17)

在求解控制方程时，采用的入口边界条件为压力梯度为0，指定入口速度；出口边界条件为压力为0，速度梯度为0；壁面处采用无滑移边界条件.由此得到伴随边界条件.

入口及壁面处

(18)

出口处

(19)

在满足伴随方程和伴随边界条件之后，式(15)剩下的部分即为灵敏度表达式

(20)

以上的推导过程采用了冻结湍流的方法，忽略了湍流黏度的变化，因此灵敏度求解结果对于层流来说是准确的，而对于湍流流动，得到的灵敏度是近似的.蜗壳优化算例结果表明，冻结湍流方法得到的灵敏度具有足够的精度，能够满足蜗壳优化的要求.

1.4 优化流程

通过改变材料的分布来改变不同单元中流体所受的阻力，流体的流速及压力分布也随之有所改变，进而可以影响目标函数的值.由此，对设计域进行优化就可以获得在一定约束条件下的最优材料分布，即蜗壳的形状.蜗壳的优化流程为

1) 对初始设计变量场进行密度过滤;

2) 通过OpenFOAM求解控制方程和伴随方程，获得速度、伴随速度等数据;

3) 计算目标函数对于设计变量的灵敏度;

4) 将目标函数和灵敏度传递给优化程序，采用MMA优化器求解后获得新的设计变量值并判断是否收敛，收敛条件为相邻2次优化计算中目标函数的相对误差不超过10-5，或优化达到最大迭代次数200次;

5) 若未达到收敛条件，则将原设计变量场更新为MMA优化器求解后获得的设计变量场，并返回第一步重新进行密度过滤及后续计算，直到达到收敛条件.

2 优化算例介绍

2.1 设计域构型

根据单级单吸清水离心泵的设计参数构建合适的蜗壳设计域，泵的质量流量Q为55.56 kg/s，扬程H为18.5 m,转速n为1 480 r/min.

蜗壳的螺旋部分是蜗壳内产生水力损失的主体部位，因此对螺旋部分单独进行拓扑优化设计.采用的二维蜗壳设计域构型如图3所示.蜗壳的基圆直径为0.27 m，底边长及总高度为0.47 m，蜗壳出口位于入口圆的正上方(蜗壳的第8断面处)，出口总高度为0.1 m，出口中心高度h0为0.05 m.约束条件中流体域的体分比上限m设置为0.45.

图3 蜗壳设计域构型

另外，在某些外部空间受限条件下，蜗壳的出口位置可能需要一定的升高，因此构建了图4中的2种蜗壳设计域进行优化设计，其出口中心高度h分别设置为0.15 m和0.25 m，体分比上限m分别设置为0.40和0.35，其他几何参数不变，以此验证该蜗壳拓扑优化设计方法的灵活性和通用性.

图4 出口中心高度升高的蜗壳设计域构型

2.2 流体流动状态

采用雷诺平均方程(RANS)下的k-ε湍流模型来描述流体的流动状态.k和ε分别代表湍流动能场和湍流耗散率场，其计算公式如下

(21)

(22)

式中:v′1，v′2，v′3分别为流体速度在x，y，z方向的脉动速度；Cμ为k-ε模型的标准参数，其值为0.09；l为湍流尺度量.

蜗壳入口速度的轴面分量为1.2 m/s，圆周分量为21 m/s，计算得到k和ε的初始值分别为1.098 m2/s2和18.905 m2/s3.在壁面上，k和ε分别采用kqRwallFunction和epsilonWallFunction边界条件.

3 优化结果

3.1 灵敏度分析结果

采用差分法验证灵敏度.随机选取设计域的3个单元分别进行灵敏度的验证，验证结果如表1所示.由于求解灵敏度时采用了冻结湍流假设，因此湍流情况下的灵敏度求解结果与差分法结果有一定偏差.

表1 灵敏度验证结果

3.2 优化结果

对蜗壳算例进行拓扑优化计算，最终得到蜗壳出口中心高度h分别为h0，3h0和5h0的设计域的优化结果如图5所示.

图5 不同出口中心高度的蜗壳优化结果

3.3 蜗壳传统设计结果

在同样的设计参数下，采用传统的一元理论蜗壳设计方法[15]进行蜗壳设计.蜗壳的基圆直径与图4中蜗壳优化设计域的参数相同，为0.27 m.

根据离心泵设计参数确定泵的比转数ns

(23)

式中：Qv为泵的体积流量.

计算得泵的比转数ns=143，根据比转数与隔舌安放角的关系[15]确定蜗壳隔舌安放角φ0=26.9°，计算得蜗壳的隔舌螺旋角α0=3.27°.

蜗壳断面形状采用梯形，根据速度系数法确定蜗壳断面1—8的面积.

蜗壳扩散管的进口近似为蜗壳第8断面，根据计算得到扩散管排出口径Dd=0.13m.根据扩散管进口当量直径和排出口径确定扩散管高度L=0.23 m，并计算得到扩散角θ=7.8°.

根据以上确定的蜗壳各部分结构参数，得到二维蜗壳构型如图6所示.

图6 传统设计方法下的蜗壳构型

将传统设计方法下的蜗壳第8断面以后的扩散管部分去掉.根据蜗壳拓扑优化结果进行精细建模，两者采用与蜗壳优化设计相同的边界条件、湍流模型等进行流场及目标函数值计算.2种蜗壳的构型及速度场对比如下.

图7 2种设计方法下的蜗壳构型对比

由以上对比可以看出，2种设计方法下的蜗壳构型具有相似性，且通过拓扑优化得到的蜗壳流动更均匀，因此目标函数值与通过传统设计方法得到的蜗壳相比较小，表明在相同的设计条件下，蜗壳拓扑优化设计方法可以替代传统的蜗壳设计方法，获得较为合理的蜗壳构型.对于出口中心高度分别为3h0和5h0的蜗壳设计域，传统的一元理论蜗壳设计方法没有明确指出蜗壳断面之间应如何过渡，无法进行此类蜗壳设计；而拓扑优化方法能够很好地适应设计域构型的改变，快速得到符合设计要求的蜗壳构型.

图8 2种设计方法下的蜗壳速度场对比

4 结 论

1) 将基于连续伴随的拓扑优化方法应用到了蜗壳的设计中，以总压损失为目标函数，推导了所对应的伴随方程和伴随边界条件，在OpenFOAM中实现了用连续伴随法求解目标函数对于设计变量的灵敏度，在一定的体积约束下，得到了二维蜗壳的优化构型.从结构来看，用基于连续伴随法的流体拓扑优化方法得出的蜗壳结构与用传统设计方法设计出的蜗壳结构十分接近；且拓扑优化方法下的蜗壳流场内的总压损失更小，由此可证明此蜗壳设计方法的可行性；同时也对不同出口中心高度的蜗壳进行了优化，最终快速得到了对应的二维蜗壳构型，验证了此方法的灵活性与通用性.

2) 目前，虽然蜗壳的传统设计方法已经十分成熟，但是对于一些在外部空间受限条件下工作的蜗壳，传统设计方法设计经验不足.而蜗壳的流体拓扑优化方法由于其不依赖于设计经验，且能够快速得到指定设计条件下的优化设计构型，因此在未来的工作中，可以被应用于外部安放空间受限以及其他复杂设计条件下的蜗壳设计中.