基于鱼形颗粒群追踪的贯流泵鱼类撞击损伤特性研究

2023-11-23张德胜史科航施卫东

农业机械学报 2023年10期

张德胜史科航潘强施卫东

(1.江苏大学流体机械工程技术研究中心, 镇江 212013; 2.南通大学机械工程学院, 南通 226019)

0 引言

泵站在防洪、发电、灌溉、跨流域调水等方面有着不可替代的作用,然而也会对生态环境造成负面影响[1]。鱼类在通过泵站时,与泵叶片碰撞造成高比例损伤或死亡,影响生物多样性并导致局部生态环境污染;水库大坝截断江河,阻隔了鱼类洄游通道,影响鱼类的产卵繁殖[2]。近年来,国家越来越重视生态环境的保护问题,《中华人民共和国长江保护法》和《十四五水利科技创新规划》明确指出保护水生生物及其洄游通道、研发水利工程过鱼设施关键技术的重要性。由此,研究鱼体撞击损伤规律并提高泵站过鱼的存活率对生态环境具有重要意义。

鱼类通过水力机械时造成的损伤和死亡存在多种因素。研究人员通过大量的实验研究发现,造成鱼类损伤和死亡的主要因素是压力波动、流体剪切力和机械损伤[3-5]。文献[6]通过实验发现不同鱼种承受压降的能力,有鱼鳔的鱼类存活率较高。文献[7]通过实验得到了鱼类可承受的剪切速率阈值为500 s-1,且受鱼种和鱼体朝向的影响。文献[8-10]通过CFD(计算流体力学)发现压力波动及剪切速率仅在叶片前缘、转轮室壁面等局部区域会对鱼体造成损伤,证明叶片撞击是鱼类损伤死亡的最主要因素。文献[11-12]进行大量的鱼类与叶片的撞击实验,发现鱼类损伤会随着叶片加厚以及撞击速度降低而减小。并建立了撞击概率与鱼体运动、鱼体长度和叶片关系的模型。文献[13-15]进行了大量的数值计算和活鱼实验,结果表明鱼体与半圆形和加厚的叶片前缘碰撞时,可以使鱼类先发生大的形变从而减弱撞击力。随后文献[16]通过鱼与叶片的撞击实验,发现叶片前缘厚度、叶片前缘倾角、撞击速度以及鱼受到的撞击部位都是鱼受到撞击后影响存活率的重要因素。

通过实验研究鱼类撞击损伤的操作复杂,实验变量不易控制且成本高昂,通过数值模拟的方法来研究鱼类经过水力机械的撞击损伤与运动特性逐渐普遍。文献[17-18]通过CFD技术与水轮机流场耦合模拟实验,发现可以通过控制鱼体进入流场入口的位置,来控制鱼体通过水轮机叶片位置,鱼体朝向是影响叶片撞击模型的最主要因素,使用随机的鱼体朝向来计算撞击概率可以提高结果的准确性。文献[19]首次通过DEM(离散元法)软件将鱼体简化为柱状颗粒,来模拟鱼体在潮汐轮机中的运动,统计鱼体的撞击概率,并引入回避率修正撞击概率结果。文献[20]通过沉浸边界和流固耦合方法研究鱼类运动,结果表明鱼体撞击损伤会随泵站流量运行的增大而增大。文献[21]通过流体力学计算与Actran软件结合研究发现,轴流泵低流量运行易导致鱼的听觉系统受损。文献[22]通过CFD模拟方法研究混流式水轮机流道对鱼类的损伤,结果表明压力损伤和剪切损伤概率与流量呈正相关,压力损伤概率是主要原因。

本文基于CFD-DEM耦合模拟的方法,通过修改编译传统的耦合接口代码,将鱼体质点当作中心点,在其周围选取多个流场点,采用矢量叠加的方法进行曳力计算。研究贯流泵中鱼体与叶片及壁面撞击后的运动行为及损伤特性,分析鱼类撞击死亡的影响因素,以期为鱼类友好型水力机械优化设计提供参考依据。

1 数值模拟方法

1.1 流体控制方程

在模拟过程遵循流场内质量、动量和能量守恒,不考虑能量传递和耗散所涉及的温度变化。流体视为不可压缩连续介质,采用RANS方程求解。流体的控制方程表示为

(1)

(2)

式中ρf——流体密度,kg/m3

t——时间,s

p——压力,Pa

g——重力加速度,m/s2

μf——流体动力粘度,Pa·s

μt——流体湍流粘度,Pa·s

F——其他作用力作用的合力,N

1.2 鱼形颗粒运动方程

鱼体在流场中不会自主运动,只受到流场的作用力。根据牛顿第二运动定律,鱼形颗粒群求解方程式表示为

(3)

(4)

式中mp——颗粒质量,kg

up——鱼体颗粒速度,m/s

Fd——流场对鱼体颗粒作用的曳力,N

Fg——粒子重力和浮力的总和,N

Fc——鱼类颗粒群间撞击产生的接触力,N

Ip——转动惯量,kg·m2

ωp——颗粒角速度,rad/s

Tf——流体力转矩,N·m

Tt——切向力转矩,N·m

Tr——滚动摩擦力转矩,N·m

Fd根据文献[23]提出的非球形颗粒的曳力模型来计算。但该公式只适用于颗粒体积小于网格体积的情况,以颗粒质心点所在网格的流场速度来代替颗粒受到的流场速度。由于本文颗粒模型尺寸远大于网格尺寸,且在流场中的运动方向随机,此时只选用颗粒质心所在网格的流场速度不能代替颗粒模型受到的曳力。本文通过修改耦合接口代码,在颗粒模型周围选取多个流体速度。如图1所示,以颗粒质心为中心,沿x、y、z正负轴以0.5倍和0.25倍鱼体的长度取点,加上颗粒质心处点共取13个,用这些点的流场速度通过矢量叠加取平均值的方法代入非球形曳力模型计算得出曳力。

图1 优化曳力计算取点示意图Fig.1 Schematic of optimal drag calculation points

(5)

(6)

(7)

b1=exp(2.328 8-6.458 1φ+2.448 6φ2)

(8)

b2=0.096 4+0.556 5φ

(9)

b3=exp(4.905-13.894 4φ+18.422 2φ2-
10.259 9φ3)

(10)

b4=exp(1.468 1+12.258 4φ-20.732 2φ2+
15.885 5φ3)

(11)

式中up——鱼体颗粒在流场中的速度,m/s

ufi——鱼体颗粒周围选取的多个流场速度,m/s

n——计算曳力时选取的流体速度点数,取13

dp——颗粒直径,m

φ——颗粒球形度,球形颗粒表面积与同体积非球形颗粒表面积之比

Fg为鱼体颗粒在流场中受到的重力和浮力总和,计算公式为

(12)

式中ρp——鱼体颗粒密度,kg/m3

Fc为鱼类颗粒群间撞击产生的接触力,包括法向分量Fcn和切向分量Fct。由文献[24]提出的软球接触模型求解,即

Fc=Fcn+Fct

(13)

其中

Fcn=-knδnn-γn(urn)n

(14)

Fct=-ktδtt-γt[(urt)t+(ωpiri-ωpjr)]

(15)

式中n——法向单位向量

t——切向单位向量

r——颗粒质心到碰撞接触点的矢量

ur——碰撞中颗粒i和颗粒j的相对速度

k——颗粒弹性刚度

γ——阻尼系数

δ——碰撞对之间的位移

流体对旋转颗粒的转矩Tf使用文献[25]的表达式,旋转系数CR由文献[25-26]的直接模拟得到,计算公式为

(16)

(17)

(18)

(19)

式中ωf-p——流体和颗粒之间的相对角速度

CR——滑移-剪切升力系数

ReR——旋转雷诺数

切向力转矩Tt、滚动摩擦力转矩Tr计算公式为

Tt=rFct

(20)

(21)

Fn=-knδnn

(22)

式中μr——滚动摩擦因数

R——与鱼形颗粒等体积的小球半径

1.3 流场-鱼体耦合

为了模拟鱼类在流场中的运动过程,通过欧拉-拉格朗日模型进行双向耦合,使用RANS方法处理流体相,使用DEM方法来描述颗粒的运动,捕捉每个颗粒运动和碰撞等信息。计算过程可归纳为:Fluent先进行单个时间步的计算,EDEM通过耦合接口,获取Fluent中的流场信息,计算颗粒所受的流体力、转矩、碰撞力等,更新颗粒运动状态,再将计算得到的颗粒位置、反作用力等通过耦合接口传递给Fluent,Fluent根据这些信息继续进行时间步计算,从而形成DEM颗粒与Fluent流场信息的相互传递。

模拟过程的具体设置如下:采用Fluent软件求解模型内的不可压三维定常流动,EDEM软件模拟鱼体颗粒在流场的运动受力,通过耦合接口实现流场、鱼体运动信息的相互传递。Fluent中采用标准k-ε模型计算,模型进口设置为速度入口,假定5%的中等湍流强度,出口设置为压力出口。模拟中将固体壁面设置为光滑、无滑移壁面。在EDEM中选择 Hertz-Mindlin(no slip) with RVD Rolling Friction 作为颗粒与壁面间以及颗粒间的相互作用计算模型。耦合设置中选用Euler-Lagrangian作为耦合方法,选用修改的曳力模型作为鱼体曳力计算方法,设置sample points值为100,增加模拟过程的稳定性。

2 物理模型与实验验证

2.1 几何模型建立

文献[16]开展了鱼类平板撞击的存活率实验,设计有5个不同前缘倾角的特殊叶片,通过控制该叶片以不同速度、不同叶片前缘倾角撞击被麻醉后用细线固定好头尾的鱼体,在叶片撞击后,鱼体可以在水箱中自由运动。以此来分析不同前缘倾角不同速度的叶片撞击后鱼的存活率,并做了统计分析表。本文的数值模拟以文献[16]的撞击实验为基础,构建一个长6 m、宽4 m、高1.5 m的长方体,内置90°、75°、60°、45°、30°共5个不同前缘倾角的叶片,叶片的前缘为10 cm厚的半圆形,如图2、3所示。在数值模拟中通过控制鱼体速度、撞击位置来进行鱼类撞击存活率预测,采用文献[16]撞击实验的部分结果拟合导致鱼类死亡的撞击力模拟阈值,用另一部分实验结果做撞击存活率模拟结果的验证。在撞击实验中,用细线固定鱼体来确保撞击部位是鱼身体的中间部位,本文的模拟实验通过EDEM颗粒工厂来初设鱼体位置和朝向,鱼体在流场中受液流作用可以自由移动,在撞击前由于只受流场流向的作用力,在撞击时保证了撞击位置和朝向与实验值相同。

图2 几何模型示意图Fig.2 Geometric model diagram

图3 网格模型示意图Fig.3 Grid model diagram

研究表明,叶片前缘越厚,鱼的变形量越大,越能提高鱼受到撞击后的存活率[13]。本文选取鱼体长度与叶片前缘厚度的比值L/d=2,鱼体长度为20 cm,与实验保持一致。使用EDEM软件通过组合多个球体建立鱼体外轮廓三维模型,由132个半径为0.5 cm球形颗粒捏合而成,见图4,鱼体模型的长、宽、厚度之比为8∶2∶1,为了达到鱼体悬浮状态,其密度与液流相同,均为1 000 kg/m3。鱼体与壁面的碰撞参数简化处理为橡胶和钢,设置如下:鱼体与叶片或壁面间的撞击恢复系数为0.95,静摩擦因数为0.63,滚动摩擦因数为0.02。

图4 鱼体颗粒模型(L=20 cm)Fig.4 Fish body particle model(L=20 cm)

对标实验进行3组模拟,设置鱼体颗粒与液流速度相同,分别以速度7、10、12 m/s从进口向5个前缘倾角不同的叶片运动。Fluent时间步长为 10-4s, EDEM时间步长设置为10-5s,来满足耦合过程时间步长的设置要求。在流场入口处,每秒生成5个鱼体,鱼体位置对应5个叶片倾角且不互相影响,数值模拟总时长30 s,每个叶片前缘可提取超过100次的鱼体撞击数据进行分析。

2.2 鱼体碰撞过程

图5为撞击速度10 m/s时,鱼体与叶片前缘倾角90°的碰撞过程。该图描绘了鱼体在流场中游动时,整个碰撞过程前后速度、曳力、撞击力和合力的变化。鱼体刚进入流场时,流场对鱼体的曳力和鱼体受到的总力很小,鱼体在流场中匀速运动。在1.795 s,鱼体的速度发生陡降,由10 m/s降到0.5 m/s,并且有撞击力出现,根据这些因素可以判断小鱼在1.795 s与叶片发生碰撞,此时出现的撞击力是小鱼与叶片第1次碰撞的撞击力,也是最大撞击力。

图5 鱼体撞击仿真性能参数变化曲线(10 m/s,叶片前缘倾角90°)Fig.5 Simulation performance parameter variation curves of fish body strike (10 m/s, 90° blade angle)

图6为撞击速度10 m/s时,叶片前缘倾角30°与鱼体的撞击过程。通过观察不同叶片前缘倾角与鱼体的撞击过程,可以发现随着叶片前缘倾角减小,鱼体受到的叶片撞击力显著减小。且随着叶片前缘倾角减小,鱼体受到的撞击情况变得复杂,一条鱼体会与叶片前缘发生多次碰撞。观察图6中1.659 s,鱼体在第3次与叶片前缘碰撞时,受到的叶片撞击力变小,鱼体速度数值几乎不变,受到的曳力却显著增大。在撞击后鱼体速度继续下降,达到谷值后回升,曳力在这个过程中持续减小。分析上述现象可得,该鱼体颗粒在第3次撞击后,速度矢量发生明显变化,然后在流体曳力的作用下,速度慢慢恢复。通过分析可得,叶片撞击对鱼体的运动轨迹有显著影响,在叶片前缘倾角小于45°的情况下,碰撞还会改变小鱼的速度方向,小鱼受到的碰撞力也会显著减小。

图6 鱼体撞击仿真性能参数变化曲线(10 m/s,叶片前缘倾角30°)Fig.6 Simulation performance parameter variation curves of fish body strike (10 m/s, 30° blade angle)

2.3 撞击存活率预测

首先在液体流速10 m/s下,通过EDEM后处理软件,在鱼体与不同前缘倾角的叶片撞击过程中,对鱼体受到的最大撞击力进行提取。以某个数值的撞击力为阈值,当撞击力小于该力时,判定鱼体在此次撞击模拟中存活,得出此力下的鱼体撞击存活率。图7为撞击速度10 m/s、叶片前缘倾角60°下模拟得到的100条鱼体与叶片前缘的碰撞力极值,若设置红色线为撞击力阈值,则鱼体撞击后存活率为黑色点数占总点数的百分比。

图7 鱼体碰撞力极值散点图(10 m/s,叶片前缘倾角60°)Fig.7 Strike force extremum of fish (10 m/s, 60° blade angle)

以不同的撞击力为撞击存活阈值,计算出在该力下鱼体的撞击存活率。再通过与文献[16]撞击实验中L/d=2、撞击速度10 m/s的鱼类存活率进行拟合,见表1,取与文献[16]实验结果方差最小的撞击力,作为导致鱼类死亡的模拟撞击力阈值。撞击力模拟阈值取变化范围2 000～3 500 N,分别得到对应的鱼类存活率预测值,并与实验撞击存活率进行方差计算,见图8。可以看出,方差的变化趋势存在极小值0.008,此时对应的撞击力阈值为2 446 N。

表1 Amaral 10 m/s实验结果Tab.1 Amaral 10 m/s experimental result

图8 撞击存活率预测值方差Fig.8 Strike survival prediction variance

以撞击力2 446 N为导致L/d=2鱼类死亡的阈值,对Amaral其他实验条件下的鱼类死亡率进行模拟预测,包括不同撞击速度、不同叶片前缘倾角及鱼体长度与叶片厚度之比L/d。将模拟结果与Amaral实验结果进行对比,见表2。

表2 模拟存活率与实验值对比Tab.2 Simulated survival rate compared with experimental value

对比表2可以看出:数值模拟对鱼类撞击叶片的存活率预测值与文献[16]实验结果保持较高的一致性,误差均在5%以内。实验和模拟结果均表明:减小叶片前缘倾角、减小鱼体与叶片的相对撞击速度、增大叶片前缘厚度可以减小鱼类与叶片撞击时的撞击力,从而提高鱼类通过贯流泵过流部件的撞击存活率。综上,通过CFD-DEM耦合模型模拟鱼类运动,得到导致L/d=2鱼类死亡的撞击力阈值2 446 N具有一定的合理性,为下文贯流泵中预测鱼类通过存活率提供理论支持。

3 鱼类过泵撞击损伤特性

3.1 水力模型

本节采用的模型泵是以某泵站水力模型为基础,通过减少叶片数、增加叶片前缘厚度和弯掠叶片的鱼类友好型设计,得到的生态友好型叶轮[27]。贯流泵模型如图9a所示,流体区域包括5部分:进水流道、叶片、导叶、灯泡体、出水流道。为保证流动的充分发展,设置进、出口延伸段长度均为叶轮外径的4倍,水泵具体参数见表3。叶轮通过弯掠叶片前缘设计来形成小于 90°的叶片前缘倾角,从而减小鱼类通过叶轮过流部件的撞击损伤。贯流泵的过鱼模拟采用CFD-DEM耦合计算,所用模型网格如图9b所示,考虑到模拟精度以及计算资源的合理分配,网格数为9.31×106,网格无关性验证在文献[27]中已开展。

图9 贯流泵模型与网格Fig.9 Tubular pump model and grid

在贯流泵的设计工况下对鱼体进行鱼类过泵撞击损伤模拟。以叶轮旋转0.25°作为一个时间步长,在流场入口处,设置EDEM颗粒工厂每秒生成10个鱼体,符合在Euler-Lagrangian 法中固相体积分数小于 10%的要求。鱼体颗粒从进口域均匀随机地跟随流体进入贯流泵流场运动,运动至出口域出口处自动移除。

3.2 鱼体运动轨迹

图10为两个鱼体在贯流泵过流部件的运动轨迹,红色是受到叶片前缘撞击的鱼体,蓝色是未受到叶片前缘撞击的鱼体。从图10可以看出,在进入叶轮区域之前,鱼体随液流运动,重力和浮力相互抵消,只受到流场力作用,近似一条直线。当鱼体进入叶轮之后,未发生撞击的鱼体运动轨迹基本与液流流线保持一致;而受到撞击的鱼体速度矢量发生改变,产生径向和周向运动,增加鱼体与各部件壁面发生碰撞的概率,进一步导致鱼体的损伤。

图10 过流部件鱼体运动轨迹Fig.10 Motion path of fish body through current component

图11为鱼体与贯流泵叶片前缘的撞击过程。在t1时刻鱼体随液流运动,在t2时刻鱼体与叶片前缘发生碰撞,运动轨迹发生变化,由于鱼体质心在碰撞点左侧,鱼体向叶片吸力面翻转,t4和t5时刻鱼体翻转后继续与叶片表面发生多次碰撞。

图11 叶片前缘撞击鱼体过程Fig.11 Process of impingement of leading edge of blade on fish

3.3 不同过流部件鱼类撞击受力分析

在贯流泵的设计工况下,选取3条长度10 cm的鱼体的过泵运动来分析鱼类通过不同过流部件的撞击过程、受力损伤、流场曳力、速度等变化,鱼体长度与叶片前缘厚度的比值L/d取2。3条鱼体分别与叶轮过流部件、导叶过流部件、灯泡体过流部件发生撞击,将受到撞击的鱼体分别命名为鱼a、鱼b、鱼c。

图12a为鱼a与叶轮过流部件叶片前缘的碰撞过程。图中出现两次撞击力峰值,表明鱼a与叶轮过流部件的叶片前缘发生了两次撞击。撞击过程中鱼a受到的撞击力最大值为2 533 N,已超过上文所得L/d=2鱼体撞击力死亡阈值2 446 N,因此可以判定鱼a与叶轮叶片撞击产生的撞击损伤严重导致鱼体死亡。

图12 鱼体与贯流泵撞击过程仿真性能参数变化曲线Fig.12 Simulation performance parameter variation curves of impact process between fish and tubular pump

图12b为鱼b与导叶过流部件的碰撞过程。鱼体b在导叶过流部件运动过程中,只在3.93 s出现撞击力的峰值1 168 N、速度的陡降、曳力的陡增,说明整个撞击过程中鱼b与导叶过流部件只发生了一次撞击,且撞击力较小,撞击力导致的撞击损伤不致命,鱼体b通过导叶过流部件仍保持存活状态。

图12c为鱼c与灯泡体过流部件的碰撞过程。鱼体c与灯泡体过流部件撞击过程中发生多次碰撞,轨迹变化也最为明显,但由于撞击速度较小受到的撞击力最小,鱼体损伤程度最小。

通过比较3条鱼在不同过流部件的撞击力可以发现,通过叶轮过流部件的鱼体受到的撞击力最大,撞击损伤最严重,通过导叶过流部件鱼体运动受到的撞击力次数最少,通过灯泡体过流部件鱼体受到撞击力最小,鱼体损伤最小。

3.4 不同过流部件鱼类撞击损伤

以相同的鱼体模型,在EDEM颗粒工厂设置等比例缩放来控制鱼体长度。在鱼体长度10 cm、L/d=2条件下进行过泵损伤模拟,并统计所有鱼体的受力信息。通过上文模拟得出在L/d=2条件下鱼体的撞击死亡阈值为2 446 N,以该阈值为判定标准,统计50条鱼在不同的过流部件中发生撞击的条数及占比,以及撞击产生后受力小于2 446 N的占比,见表4。得出模拟中L/d=2的鱼体通过贯流泵叶轮过流部件撞击概率为18%,撞击存活率为66.7%,综合鱼类存活率为94%,鱼类通过其他过流部件的碰撞死亡率近乎为0。从撞击力小于2 446 N的占比来看,叶轮造成的鱼类损伤最严重,导叶和灯泡体造成的鱼类碰撞损伤较低。

表4 鱼体通过贯流泵的碰撞比例 (50条鱼)Tab.4 Strike ratio of fish through tubular pump (50 fish)

3.5 不同L/d鱼体撞击损伤

为了分析在不同L/d条件下的鱼体撞击受力,通过不同鱼体长度来改变L/d的比值,图13和图14分别为50条鱼体在L/d=4和L/d=8的条件下,通过叶轮、导叶和灯泡体3个过流部件的撞击位置、撞击力极值统计图。通过比较不同过流部件的撞击力,发现鱼体在叶轮区域受到的平均撞击力最大,在灯泡体区域受到的平均撞击力最小,叶轮撞击对鱼类存活的威胁最大。鱼体通过叶轮区域时,发生碰撞的位置包括鱼体头部、腹部或尾部与叶片前缘、叶片表面、叶轮轮毂撞击等,结果表明:鱼体与叶片前缘撞击时,鱼体受到的撞击损伤最高;L/d=8的鱼体受到叶片前缘的撞击力高于L/d=4的鱼体,说明L/d的变大会导致更高的鱼体撞击损伤。

图13 L/d=4的鱼体通过不同过流部件撞击力Fig.13 Strike force of fish body with L/d=4 through different flow parts

图14 L/d=8的鱼体通过不同过流部件撞击力Fig.14 Strike force of fish body with L/d=8 through different flow parts

统计L/d=4和L/d=8的鱼体通过不同过流部件发生的撞击次数,见表5。可以发现L/d=8的鱼体在所有过流部件中的撞击次数均高于L/d=4的鱼体,L/d的取值与撞击概率显著相关[28]。在叶轮过流部件中,L/d=8鱼体的碰撞力超过2 446 N的碰撞次数占比为40.3%,远高于L/d=4鱼体的占比16.1%,L/d变大,鱼体碰撞受力超过2 446 N的占比变大。此外结合图13、14发现,L/d=8的鱼体在叶轮叶片前缘处发生撞击次数接近L/d=4的鱼体的两倍,叶片前缘撞击次数与L/d的比值成正比[10]。