四川工业学院 学报
Journal of Sich uan Uni v ers ity o f Sc ience and Technology
文章编号:1000 5722(2000)02 0079 06
收到日期:1994 11 18 国家博士点基金资助项目
作者简介:刘小兵(1965-),男,四川阆中人,四川工业学院能源与环境工程系系主任、教授、工学博士,主要从事水力机械的汽蚀与泥沙磨损、多相流动理论和数值模拟研究。
水力机械泥沙磨损的数值模拟
刘小兵
(四川工业学院能源与环境工程系 四川成都610039)
摘 要: 作者对低浓度沙水(平均固体容积浓度C v <1%,或含沙量低于20kg/m 3)中运行的水力机械的泥沙磨损进行了数值模拟研究。利用Euler ian-lagrang ian 混合湍流模型(这个模型是基于两相中的颗粒源项数值技术,即液相采用Euler ian 描述,而颗粒相采用Lagr angian 描述)对水力机械中湍流速度、压力分布以及颗粒轨迹进行了数值计算。同时给出了水力机械中的颗粒反弹模型和磨损模型。通过对一水轮机导叶磨损率的计算,其结果与实际较为吻合。
关键词: 数值模拟;泥沙磨损;水力机械 中图分类号:TV 131;T K72
文献标识码:B
我国河流的特点之一是含沙量比较大,就长江来说,年平均输沙量已达5 14亿t,三峡库区年平均含沙量1 17kg/m 3,最大达10 5kg/m 3。多泥沙河流上水电站水轮机的磨损问题是客观存在的自然现象,为国内外水轮机界所重视,尤其是我国三峡巨型电站现在已经上马。尽管国外的许多电站含沙量还不到0 1kg/m 3,但其水轮机过流部件磨损仍很严重。因此,预测和减轻乃至解决含沙水流中水轮机的磨损问题是非常必要的,也是当今的重大课题。从影响水轮机泥沙磨损的主要因素可知,掌握沙粒在过流通道内的运动规律对研究水轮机的泥沙破坏规律起着重要的作用。由于场中的沙粒运动受到各种力、湍流以及空化的影响,使得含沙水流的流动规律非常复杂,因此,目前对各种材质的抗磨性能、磨损机理等基础性试验研究开展较多。近年来,随着计算机计算容量增大、计算速度提高,沙水对过流部件磨损的数值模拟(数值实验)迅速发展。
针对水轮机中含沙水流的浓度并不很高,本文建立了低浓度固液两相流Eulerian -Lagrangian 混合湍流模型,平均固体体积浓度C v <1%,即对于沙水,含沙量低于20kg /m 3。该模型的实质就是在Eulerian 坐标系下描述流体相的运动,在Largrangian 坐标系下处理颗粒相的运动,两个坐标系下的运动方程由两相间作用的颗粒源项来耦合,来研究颗粒群和流体相之间的较大滑移,给出了颗粒-壁面碰撞模型和为韧性金
属材质的水轮机过流部件的磨损模型。利用这些模型可数值模拟水轮机中含沙水流的流动,沙粒在水轮机过流通道内的运动轨迹及水轮机过流部件的磨损率。本文模拟了长江沙水对ZZ560-LH 型水轮机导叶部件的磨损率,并给出了其严重磨损区域。对数值模拟结果同实际进行比较,结果较一致。文中假设沙粒为球形,不考虑沙粒撞击壁面发生的破碎现象、沙粒间的碰撞及空化的作用。
1 运动方程
1 1 沙粒群运动方程
沙粒的运动在Lagrangian 坐标系中考虑,沙粒群按初始尺寸分组,各组沙粒沿其自身轨道运动。由于沙粒在场中受到阻力、压强梯度力、虚拟质量力、Basset 力、Saffman 升力、Magus 升力和重力等力的作用,因此沙粒群将沿轨道发生速度变化。在沙粒沿其自身轨道运动时,对水流造成了分布于整个体积中的物质源、动量源。
在Lagrangian 坐标系中考虑沙粒群沿轨道运动时,已由我们[1,2]修正和发展了前人所建的单颗粒的运动方程,其所建方程较全面地考虑了颗粒在场中的受力,适用于任意场中的颗粒运动,其方程表达式为:
d V p i d t =1K mi +S 34d C D
|V f -V
p |(V f i -V p i )
+
3
2d K Bi
f
t-d V fi d-d V pi d d t-
+6
d K s| f
!V f j
!x i|
1
2(V fj-V pj)+
3
4
C M!!i
?(V f i-V pi)+(1+K m)d V f i
d t
- f?2V f i -(1-S)g i](1)式中,K m(#0 5),K B(#6 0);C M(#1 0)和K S(# 1 615)分别表示虚拟质量力、Basset力、Magus升力和Saffman升力系数;d为颗粒直径;S为颗粒密度?p与流体密度?f的比值,对于沙水,一般有?p= 2650kg/m3;?f=1000kg/m3;S=2 65;g i为重力加速度分量,x i为坐标分量, f为流体运动粘性系数;!!i =!p-0 5??V f i,!!i?(V f i-V p i)表示!!?(V f-V p)的i分量,粘性阻力系数C D与颗粒雷诺数Re p有关,其表达式有:
C D=24/Re p
(24/Re p)(1+0 15Re0 687
p
)
0 44
Re p=|V f-V p|d/ f<1
1#Re p#1000
Re p>1000
脚标:f和p分别表示流体和颗粒,i和j为张量坐标。
方程(1)中Saffman升力项仅在大速度梯度场中(一般在边界层流动中)比较重要,如果不计颗粒的自由旋转速度!p,Magus升力则比Saffman升力小一个数量级,可忽略,粘性影响项- f?2V f i较小,也可忽略。
1 2 水流运动方程
固液两相湍流流动的三维平均微分方程在笛卡儿坐标系(x,y,z)中的通用形式可表示成:
!
!x(V f x#)+
!
!y(V f y#)+
!
!z(V f z#)=
!
!x?#
!#
!x+
!
!y?#
!#
!y+
!
!z?#
!#
!z+S#+S
p
#
(2a)对于不可压缩的轴对称流动,其时间平均微分方程在旋转柱坐标系(r,%,z)中的通用形式可简写成: !
!z(V f z#)+
1
r
!
!r(rV f r#)=
!
!z?#
!#
!z+
1
r
!
!r r?#
!#
!r+S#+S
p
#
(2b)式中,#为流动参数代表值;S#为源项;S p#为附加颗粒源项,表示由于流体一颗粒作用,颗粒相使流体相引起的动量变化;?#为湍流扩散值。它们在方程(2b)中的值详见表1。
表1 方程(2b)中的湍流扩散,各源项值及K &湍流模型中的实验系数值方程#?#S#
连续方程100
V fr e ff-1
?
f
!p
!r+
!p
!z r e ff
!V fz
!r+
1
r
!
!r+r e ff
!V fr
!r+
1
r(V f%+!r)
2-2
r2
eff V fr
动量方程V fe e ff-1
r V f r(V f%+2!r)-
1
r2
V f%
!
!r(r eff)
V fz e ff-1?
f !p
!z+
!
!z e ff
!V fz
!z+
1
r
!
!r r eff
!V fr
!z
K-方程K eff/?K G-&
&-方程& e ff/?&&(C1G-C2&)/K
G= e ff2!V
fz
!z
2
+
!V
fr
!r
2
+
V fr
r
2
+
!V
f%
!z
2
+
!V
fz
!r+
!V
fr
!z
2
+
!V
f%
!r-
V f%
r
2
e ff= f+C(?
f K2/&, C(#0 09, C1#1 44, C2#1 92, ?k#0 9, ? 3
1 3 附加颗粒源项S p#的计算
由于颗粒相的存在和运动会影响到流体相的流场,而流体相的流动变化也必然影响颗粒相的运动及其运动轨道,因此流体相的流动和颗粒相的运动计算必须相互耦合,才能真实地描述出流体相和颗粒相的运动情况,流体相流动的Eulerian方程式和颗粒相Lagrang ian方程式的相互耦合是通过颗粒相在每个计算单元中的源项S p#来实现的。
用Lagrang ian方法求解沙粒运动时,每一条沙粒轨迹代表一群直径相同的沙粒的运动情况。如果用一组离散为直径不同的沙粒(k)来表示沙粒尺寸的连续分布,定义f l o为起始第l站沙粒质量份数,f k0为起始处直径为d k的沙粒的质量份数,
m p0为沙粒的总流入率,则沿第l条轨迹直径为d k的运动沙粒质量流入率为
m p k=
m p0f k0f l0
而沿第l条轨迹运动,初始直径为d k的这一群沙粒的数目流量率为:
80四川工业学院学报 2000年
N pk = m p kl /m pk l =6 m p 0f k 0f l 0/( ?p k d 3
k )
沙粒在沿轨迹l 穿过第n 个计算单元时,引起的沙粒质量源项为:
M p kl = N p kl [(m pkl )out -(m p kl )in
脚标:out 和in 分别表示出去和进入计算单元时刻,而第n 个计算单元内的总的沙粒质量源项为:
(S p 1)n =
?
k
?l M
pkl n /(?f )n )=
?k
?l
N
pkl [(?p k d 3
k )out
-(?pk d 3
k )in ]n/(6?f )n )
式中,?
k
表示对通过第n 个计算单元的各尺寸组k 求
和,
?
l
则表示对通过该单元的同一尺寸组不同初始
值或轨迹l 求和,)n 为计算单元体积。
颗粒的动量源项类似地可写成:(S p vf i )n
=
?k
?l
N
p kl [(V i pkl m p kl )out
-(V i pkl m pk l )in ]n /
(?f )n )
=
?k
?l
N
p kl [(V i p kl ?p k d 3
k )out -(V i p kl ?pk d 3
k )in ]
n /
(6?f )n )式中,上脚标i 表示速度的坐标分量,湍动能源项S p k
和
其耗散率S p &可取为零。
2 沙粒浓度场的模拟
用Lagrang ian 方法来描述沙粒运动时,可求得沙粒的轨迹和在轨迹上任意点处的速度。但是在工程实际问题的应用中,往往需要知道沙粒的浓度场和速度场的分布情况,这不单是为了能把计算结果易于和实验结果相比较,而且对于分析讨论沙粒的运动和提出改进流场情况达到减小磨损率、提高运行可靠性都具有指导意义,这里引用平均体积法来求解浓度场。
在每一个计算单元内,沙粒总的数密度可由下列关系式确定:
( N )n =
?k
?l
N
kl ?t k l n /
)n 式中:?t k l =t out -t i n 是沙粒在单元内的停留时间。
则该单元内的沙粒体积浓度为:(C v )n =
?k
?l
N
k l ?t kl )p kl n /
)n
=( /6)
?k
?
l
N kl ?t kl d 3
k n /
)n
3 磨损模型
过流部件表面被固体颗粒撞击而遭磨损的问题是
很普遍的,因而目前存在很多磨损模型,针对文中的固体颗粒为长江泥沙(多为石英砂,约占65%),水轮机金属部件又多采用韧性金属材质,因此我们选择由Grant 和Tabakoff [3]所给出的单位质量的粒子碰撞固壁所产生的质量磨损率 E 的经验表达式:
E =K 1(V p 1cos +1)m 1(1 0-R 2T )f (+1)
+K 2(V p 1sin +1)m 2
式中,
R T =1-K 3V p 1sin +1
f (+1)=
[1+K 4sin (90%&+1/+0)]
1
+1#2+0
+1#2+0V p 1和+1分别为沙粒撞击固壁的速度和入射角;m 1和
m 2为速度指数;+0为最大磨损的颗粒入射角(一般在20%~30%之间),取决于材质特性;K 1,K 2,K 3和K 4分别为取决于材质特性的经验系数。针对本文的计算对象,取m 1=1,m 2=4,+0=25%,K 1=1 505?10-6[(s/m)2],K 2=
5 0?10-12[(s/m )4],K 3=
0 0016(s/m ),K 4=0 296077。
4 数值计算方法
4 1 边界条件和初始条件
边界条件和颗粒运动的初始条件(颗粒的初始位置,初始运动特性及初始浓度分布)给定得正确与否,将直接影响最后求解结果与实际情况的相符性,因此
它们的给定是很重要的。
流体(水流)相:进口边界上,给定流动的速度分布,假定流动已充分发展为湍流,则可给定湍动能K 和能耗率&为:
K |
in
=0 005(V 2f i )in &|in =3K 3/2
in /L
式中,L 为一特征长度,如对于蜗壳,则为蜗壳进口管
径,对于导水机构,则为导叶高度,对转轮,则为转轮进口边长度。
在固壁上,流体的切向和法向速度分量均为零,即V f i |
W
=0。输运方程中的湍动能K 可根据边壁应力产
生K 的表达式,临近壁区的&值根据湍流长度尺寸与来自壁面距离成线性变化来确定[4]。出口边界条件根据具体情况而定。
颗粒(沙粒)相:颗粒的初始条件,即初始速度的给定可以根据下表达式来确定:
V p 0=e 0V f o
式中,e 0为初始颗粒跟随流体的跟随系数,反映了颗粒的跟随程度。系数e 0主要与颗粒尺寸和颗粒密度有关,大颗粒和重颗粒的跟随系数e 0较小,小颗粒较大,
81
第19卷第2期 刘小兵:水力机械泥沙磨损的数值模拟
一般取为1。一般情况下,可认为在进口处颗粒的速度均匀分布,当然也可以设颗粒速度按一定规律分布,并且不同尺寸组的颗粒具有不同的初始速度,这将更近于实际情况。为了获得一个光滑的浓度曲线,应该给定足够多的颗粒位置的计算站。在进口处,通常用离散直
径来表示颗粒尺寸的连续分布。离散直径取几档应该根据具体情况确定,一般取3~5个尺寸段,以充分描述分散颗粒群的运动规律。
沙粒与固壁的作用采用颗粒-壁面碰撞模型,即颗粒撞击壁面的入射角+1可由下式来计算:
+1=tg -1(V p n 1/V pt 1)
当颗粒撞击壁后,反弹速度的切向和法向分量可根据Tabakoff 和H amed [5]中石英砂与不锈钢板碰撞的反弹关系式:
V pn 2/V p n 1=1 0-0 4159+1-0 4994+21
+0 292+31
V pt 2/V p t 1=1 0-2 12+1+3 0775+21-1 1+3
1
颗粒撞击固壁后的速度V p 2和反弹角+2由下式来确定:
V p 2=
V 2pn 2+V 2p t 2 +2=tg -1
(V p n 2/V p t 2)
式中,V pt 和V p n 分别为颗粒撞击点的切向和法向速度,脚标1和2分别表示碰撞前后。4 2 数值计算方法
首先应该求解出流体相的流动,本文将三维导叶图1 计算网格示意图
流场分解成两个二维流动,即平面叶栅流动(沿导叶高度方向均匀取8个)和轴面流动,并使用SIMPLE 方法[6]对流体相的两个二维流动进行相互迭代求解,计算网络的划分见图1(平面叶栅流场有159?24个结
点,轴面流场有167?32个结点)。在计算的第一步,可不考虑颗粒源项,即计算出无沙粒时的水流场,以此作为两相流计算的初始场,在所得初始场下求解沙粒运动方程。在得到了流体相的流动特性后,通过数值求解颗粒运动方程(1),求出颗粒速度、轨迹、颗粒群源项,如果颗粒与壁面相撞,则应使用颗粒-边壁碰撞模型,以计算出颗粒碰撞固壁的入射角,反弹角和反弹速度,然后把所求出的颗粒源项代入流体相方程再进行求解,修正水流场,以所求得修正后的水流场再求解颗粒相方程,如此迭代计算下去,直至流体相与颗粒相分别收敛为止。最后利用这些流动参数和磨损模型求出水轮机导叶部件的磨损率,计算的流程框图如图2所示。
颗粒运动方程(1)的求解一般采用四阶Runge Kutta 积分法,也可用差分等其它方法,在求得沙粒的运动速度后,就可利用微分方程
d x p i
d t
=V p i 求出沙粒的运动轨迹。上述方程也可采用四阶Runge Kutta 方法,不过一般采下列近似方法求解:
x p i (K )=x p i (K -1)+0 5[V p i (K -
1)
+V pi (K )]?t
式中,?t 为时间步长;K 为计算点,V p i (K )为时刻K ?t 的沙粒速度。
5 数值计算结果和分析
作者选用的沙样取自葛洲坝电站过机泥沙,粒度
82四川工业学院学报 2000年
分析见表2。选用机组为葛洲坝ZZ560-LH 机型,计算模型的转轮直径均为420mm,导叶材料为ZG25。导叶布置如图3。实验情况表明设计工况下的最大磨损率比非设计工况下的小得多,因此这里主要针对非设计工况进行计算。
表2 泥沙粒度分析
粒径(mm )百分比(%)
1 0~7 50 3~1 08 00 15~0 310 20 075~0 1538 5~0 075
35
8
图3 导叶布置图
图4示出了直径d=0 1mm 的沙粒在大小开度下a 大开度,冲角,=10%,???沙粒反弹后的运动轨迹b 小开度,冲角,=-8%,?
?
?沙粒反弹后的运动轨迹
图4 沙粒(d=0.1mm)在导叶水平面流场中的运动轨迹
(换算后的冲角,分别约为10%,-8%),在导叶高度中部水平面流场中的运动轨迹(由于周期性边界,仅示出了单个导叶情况)。由图4可知,在大开度下导叶表面
的冲击区主要在导叶头部、导叶表面的内侧(面向转轮侧)和导叶表面外侧的尾部(由于回流漩涡卷入的沙粒的冲击)。在小开度下,导叶表面的冲击区主要在导叶头部、导叶表面外侧和导叶表面内侧尾部(由于回流漩涡卷入的沙粒的冲击)。不过由于沙粒撞击导叶表面反弹后,可能再次冲击表面,所以大小开度下,导叶表面各处都有可能遭到冲击,越大的沙粒再次冲击表面的可能性越大。越小的沙粒越易被导叶尾部回流漩涡卷入,并带有很大的旋转速度冲击导叶表面。计算结果还表明,粒径小于0 05mm 的沙粒轨迹基本与水流流线重合,在内侧运行的粒径较大的沙粒,其轨迹偏离流线向导叶表面靠近,而外侧附近的较大沙粒偏离流线远离导叶表面,粒径越大偏离越大。
图5示出了直径d=0 1mm 的沙粒在导叶轴平面流场中的运动轨迹(大小开度下的运动轨迹相差较小)。从图5可知,沙粒轨迹向导叶下端面偏离,尤其是导叶尾部。沙粒向下端面的密集造成越往下端面浓度越大,在浓度曲线图6中可以看出。计算结果还表明沙粒轨迹向上端面偏离水流流线,沙粒越大的偏离
越大。
???水流流线
沙粒轨迹
图5 沙粒d=0.1mm 在导叶轴平面流场中的运动轨迹
图6示出了平均沙粒浓度C v =0 5%,含沙量为10kg/m 3的长江含沙水流在导叶场中的浓度变化。
图6a 给出了从导叶内侧到相邻导叶外侧沿D =0 95D 0(D 0为导叶分布圆直径)圆断面上平均浓度的变化。可以看出,在大开度下,从内侧,浓度逐步减小,但到外侧稍有回升,在小开度下从内侧浓度稍下降后逐渐增高。浓度在大开度下外侧稍有回升和小开度下内侧较小下降,主要是由于导叶尾部的漩涡卷入大量沙粒所致。图6b 给出了沿D=0 95D 0圆周断面,在导叶高度方向上平均浓度的变化。可以看出,在大开度下,浓度浓渐增高,小开度的变化和大开度的基本相同,但在上端面附近比大开度的稍大,而在下端面附近稍小。
83
第19卷第2期 刘小兵:水力机械泥沙磨损的数值模拟
a 浓度沿圆周断面上的变化
b 浓度沿导叶垂直断面上的变化S 为导叶内侧计起的圆弧长,Z 0为导叶数;b 为离上端面的距离,b 0为导叶高度
???大开度,冲角(,=10%),??
?小开度,冲角(,=-8%)
图6 含沙水流在导叶场中的浓度变化
图7示出了平均沙粒浓度C v =0 5%,含沙量为10kg/m 3的长江沙水对导叶表面的平均磨损率,表明了最大磨损率在导叶头部,
随后下降。但在大开度的
a 大开度(,=10%)
b 小开度(,=-8%)???导叶内侧(模拟值), ??
?导叶外侧(模拟值)。
&导叶内侧(实验值),?导叶外侧(实验值)图7 含沙水流对导叶表面的平均磨损率
导叶外侧尾部和小开度下的导叶内侧尾部的磨损率曲线回升,其原因主要是漩涡卷入大量沙粒冲击导叶表面。大开度下,导叶内侧比外侧磨损严重,小开度则相反。计算结果还表明随着沙水浓度增高,磨损率明显
上升。偏离设计工况越远,磨损率越大,尤其是导叶头部的冲击磨损和导叶尾部的漩涡卷入大量细沙的磨损。
数值模拟结果比实验结果有所偏小,我们认为这主要是由于模型中所给定的假设所致,不过其计算结果已足以反映出过流部件的磨损率。
6 结论
针对含水沙流中运动的水轮机,建立了低浓度固
液两相流Eulerian-Lagrangian 混合湍流模型,给出了颗粒-壁面碰撞模型和为韧性金属材质的水轮机过流部件的磨损模型。用这些模型模拟了ZZ560-LH 型
水轮机导叶流场中的长江沙水流动,浓度分布及对导叶部件的磨损率,并给出了其严重磨损区域。数值模拟结果与实际结果比较一致。算例的模拟结果已表明了本文的模型和计算方法可数值模拟水轮机中含沙水流的流动、沙粒在水轮机过流通道内的浓度分布、运动轨迹及水轮机过流部件的磨损率,也可用于其它固液两相流动中。
参 考 文 献
[1]刘小兵,程良骏 固体颗粒在水涡轮机械中的运动[J ] 华中理工大学学报,1994,22(1):10~16
[2]刘小兵,程良骏 空泡在任意流场中的运动研究[J ] 水动力学研究与进展A 辑,1994,(2):150~162
[3]Grant,G ,W Tabakoff Erosion Rrediction in Turbomachinery Res ulting from Environmental Solid Particles [M ] J Aircraft,1975,(12):471~478
[4]Pourahmadi,F Turbbulen ce Modeling of Single and Two Phase Curved Channel Flow s[J ] Ph D Th esis,University of California,Berkeley,1982
[5]
Tabakoff,W ,A Hamed Aerodynamic Effect on Erosion in
Turbomachery[J].JSM E and ASME,paper ,1977,(70):392~401
[6]Patankar ,S V Numerical Heat Transfer and Flu id Flow [M ].Hemisphere,New York,1980
&YL &
Numerical Simulation of Silt Abrasive Erosion in Hydraulic Machinery
LIU Xiao bing
(Dept.of Energy and Enoi ronment Engineering of Sichan University of Scien ce and Technology Chengdu 610039)
Abstract:Numerical investigation into silt abrasive erosion of hydraulic machinery for low sand water concentration (mean sand concentration by volume C v <1%,or sand carrying capacity less than 20kg/m 3)is reported The turbulent flow velocity fields and pressure distributions as well as the trajectories of sand particles in hydraulic machinery flow are numerically calculated using the mixed Eulerian Lagrang ian turbulence model This turbulence model is based on the particel source in cell numerical technique for two phase flows,and,Eulerian equations are used for the description of the liquid phase dynamics,whilst the solid particles are treated within the Lagrangian framework The sand particle rebound model and the erosion wear model in hydraulic machinery flow field are also presented We calculated the erosion wear rates of a hydraulic turbine guide vane The numerical predictions show good agreement with the actual
Key words:numerical simulation;silt abrasive erosion;hydraulic machinery
84四川工业学院学报 2000年