摘要：本文采用变密度流基本方程和混合有限分析法，求解了突扩边界下的异重流潜入运动。通过对数值计算结果分析，讨论了进口密度弗汝德数Fre对潜入运动的影响，加深了对异重流潜入规律的了解。
关键词：异重流 潜入运动 数值计算 混合有限分析法
1 前言
异重流是自然界中常见的一种流动现象，如水库中的浑水异重流，河口地区出现的盐水楔，工业中遇到的温差异重流，底部突然下沉的小型沉沙池中的泥沙异重流等，在水库、船闸引航道、河港、河口、沉沙池等实际工程中，浑水异重流的淤积问题往往是必须考虑的问题[1]。
由于异重流潜入发生处的边界条件复杂，具有与一般明渠流不同的水流泥沙运动规律，且潜入处的水力泥沙因素与异重流的发生、发展及运行状况密切相关[2]，历来广为研究者们所关注。用一维的方法，再加上一些经验关系解异重流运动，无疑具有其实际意义，但异重流运动具有较强的三维性，不是一维流动理论基础所能较好解决的，并且在用一维方法求解异重流运动时，常常碰到诸如上下水层交界面位置及交界面上水体之间的动量和质量交换系数等难以确定的困难。近年来，逐步采用数值模拟技术对异重流潜入运动进行研究，如Yoon T.H。采用流函数涡量法对异重流的潜入运动进行了模拟和研究[3]；方春明、韩其为用立面二维对异重流潜入流动进行数值模拟，探讨了影响异重流潜入条件的因素[4]。但由于问题本身的复杂性，许多问题仍有待进一步研究。
|
考虑到异重流潜入运动的过程中，不仅具有初始动量，而且还受到因密度差引起的浮力作用，与静止环境中的浮力射流运动有相似之处，不同之处在于异重流潜入发生处的边界条件复杂，因此，本文将异重流的潜入运动视为一种具有复杂边界的浮力射流运动，从Navier Stokes方程出发，引入满足布辛涅斯克假定的一组变密度方程，并采用混合有限分析法对该运动进行一系列的数值模拟，通过对数值计算结果进行讨论和分析，进一步加深对该流动的认识。 2 数学模型 |
|
图1所示：密度为(ρ0+Δρ)的浑水从水平河道流入密度为ρ0的清水矩形水库中，由于密度分布不均匀，流场中就出现了浮力，在几何形状突然扩大的流道中，形成了有浮力回流的流动。图中x、z分别表示水平和竖直方向，进口浑水流速u0满足抛物线分布。
2.1 控制方程
该模型遵循以下假定：忽略泥沙与水流之间的速度差异；根据布辛涅斯克近似，仅重力项计及密度变化；认为泥沙在水平和竖直方向上扩散系数相同；仅对异重流潜入段进行模拟，此时浮力作用较大，且计算区间不太长，因此在这里只考虑异重流的对流和扩散，不考虑其沿程沉降[4]。
为使计算和讨论具有一般性，以水库初始水深H、进口平均流速u0为特征尺度，引入无量纲量X=x/H，Z=z/H，U=u/u0，W=w/u0，T=tu0/H，P*=p/ρu20
则控制方程的无量纲形式为
|
连续方程 |
|
|
|
(1) |
|
运动方程 |
|
|
|
(2) |
|
|
(3) |
|
浓度方程 |
|
|
|
(4) |
|
式中 Γ=ρ-ρ0/ρ0为浑水的有效密度，Re=U0H/v为雷诺数，Fr=U0/为弗汝德数，Sc=v/D为斯密特数，v、D分别为流体的运动粘性系数和泥沙扩散系数。 |
2.2 初始和边界条件
2.2.1 初始条件
初始时刻，水库水体静止，进口断面的浑水流速满足抛物线分布，且认为此断面的密度分布均匀；
初始压力按静水压力给出：P(X，Z)=1-Z。
2.2.2 边界条件
研究的流动区域共有进口、出口、自由表面、上游固壁及底部五个不同的边界，如图1所示。
进口：浑水流速满足抛物线分布，断面密度均匀分布，即
U=7.5(6-5Z)(Z-H-H0/H)，W=0，Γ=Δρ/ρ0;
|
出口：认为流动已充分发展，则; 壁面：给定无滑移边界条件，即U=0，W=0，(上游)或(底部)； 自由表面：对自由面作刚盖假定，并采用对称平面条件，即。 |
3 计算方法
3.1 求解区域及网格划分
该流动问题是开边界问题中的一种，当出口采用充分发展假定时，求解区域应保持足够的长度，以便在我们感兴趣的区域准确地模拟出流体的运动规律。依文献[5]，本文计算区域长度取为L=20。因水库进口附近的流动是主要的，远离进口的流动是渐变的，故在数值计算中，X方向上采用非均匀网格，进口附近的网格密，越远离进口，网格越疏；在Z方向上，除了0、1附近的网格密些（与边界条件有关），其余部分采用均匀网格。依控制方程及边界区域等特点，采用原始变量法求解为便，故采用交错网格来离散并求解整个流场和密度场[6]。计算时，动量方程和浓度方程的控制精度为10-7，连续方程的控制精度为10-4。
3.2 计算格式
采用混合有限分析法的五点格式进行离散。
在局部单元格上，将非线性方程（2）～（4）线性化，写成下列一般形式
|
|
(5) |
则各方程的系数可归纳为表1的形式。
方程(5)的混合有限分析解为
|
|
(6) |
表1 控制方程系数
Coefficients of governing equations
|
| ||||
|
φ |
E |
2A |
2B |
S |
|
| ||||
|
U |
Re |
ReU |
ReW |
|
|
| ||||
|
W |
Re |
ReU |
ReW |
-Re1/Fr2Γ |
|
| ||||
|
Γ |
ScRe |
ScReU |
ScReW |
0 |
|
| ||||
|
式中 Ci-1，j、Ci+1，j、Ci，j-1、Ci，j+1和Cp分别为混合有限分析法系数，具体表达见文献[7]。其中，(S-Eφt)n为第n次的迭代值，Sn用中心差分计算，φnt用向前差分表示，即()np=φij n+1-φnij /Δt，将(i-1，j)、(i+1，j)、(i，j-1)、(i，j+1)依此用1、2、3、4表示，则式（6）可写为 |
|
|
(7) |
对整个计算区域，由式（2）～（4）形成代数方程组，它们都是对角占优的，可用追赶法求解。
4 计算结果和分析
4.1 Fre对潜入形态的影响
异重流在潜入水库的过程中，不仅具有初始动量，而且还受到因密度差而引起的浮力作用，流动是非定常和变密度的，此时影响流动的一个重要的特征数是进口密度弗汝德数Fre，具体表达形式为
|
|
(8) |
式中u0为进口处浑水的特征流速，L为特征长度，在本文中，u0等于进口处浑水的平均流速，L等于进口水深H0，ρ为进口处浑水的密度，ρ0为水库的清水密度。Fre反映了进口水流惯性力作用与浮力作用之比：当其它条件不变时，Fre越小，说明浑水惯性力越小，浮力作用越大；反之亦然。
为研究异重流潜入形态随Fre的变化规律，本文在进口单宽流量一定的情况下，通过改变进口处浑水的密度ρ，分别计算了Fre=0.25，0.501，1.00，1.58时异重流各特征量的变化情况，并绘出了进口附近区域在T=0.8（T为无量纲化时间，T=tu0/h）时刻相应的有效密度等值线图（Γ/Γ0）和速度矢量图（图2和图3），由图可看出，Fre 对异重流潜入流动特性的影响是较大的。
|
图2是不同Fre的相对有效密度等值线图：当Fre=0.25时，负浮力占主导地位，浑水入库后不久便下沉，在进口附近形成一股潜流；当Fre=0.501，1.00，1.58时，入库浑水的初始动量作用逐渐增大，负浮力失去其主导作用，浑水继续向前水平运动，直到水库清水对其阻力与维持其运动的力平衡为止，此时因负浮力作用浑水开始潜入，该位置通常也称为潜入点。从该图还可以看出，Fre越大，浑水与水库水体的混合越充分，在水库中的扩散程度越大，潜入点距进口的距离越长，越不容易形成异重流；反之，则越容易形成异重流。图3是相应的速度矢量随Fre的变化图。该图反映了在进口下游一定区域内会形成较大的回流旋涡区，随着Fre的改变，回流旋涡区的形状、大小和中心位置相应发生改变。 |
|
当Fre较小时，浑水在进口附近区域很快吸附近岸，在水库中形成一个逆时针方向的旋涡，且该回流旋涡区范围较大，密度场的范围相应地也较大；随着Fre的增大，在进口附近将出现一顺时针方向的旋涡，同时逆时针方向的回流旋涡区范围将随之减小，密度场的范围相应地也减小；当Fre增大到一定程度后，进口附近区域仅出现顺时针方向的回流旋涡，进口附近水流的掺混作用剧烈。
|
|
|
|
图3 异重流潜入时的流场随Fre的变化(T=0.8) |
图4 Xp/H与Fre之间的关系 |
为进一步研究异重流潜入形态与Fre之间的关系，我们还利用数值计算的结果，在图4中点绘了潜入点与进口之间的无量纲距离Xp/H与Fre之间的关系。该图清楚地反映了两者之间的关系：Fre&<0.6时，Xp/H随Fre变化不大，进口附近的流动主要受负浮力作用的影响，此时的流动为卷流型；Fre&>1.5时，随着Fre增大，Xp/H迅速增大，进口附近的流动主要受浑水在进口处的初始动量所控制，此时流动为射流型。在卷流型和射流型之间存在一个过渡区，过渡区的流动不仅具有一定的初始动量，而且还受到负浮力作用，此时流动为浮力射流型。若不考虑过渡区，则Fre&<1时，流动为卷流型；Fre&>1时，流动为射流型，与文献[3]的结论相同。
转贴于4.2 异重流潜入的发展过程
异重流潜入水库的过程中，流动是非定常的，流场和密度场均随时间变化。为研究异重流潜入随时间的变化规律，本文根据数值计算结果，绘制了Fre=1.58的情况下，进口附近区域分别在T=0.8，1.6，2.4，3.2时刻的速度场和密度场，如图5和图6所示。
|
图5是不同时刻下的速度矢量图：T=0.8时，因初始动量的作用，在进口附近将产生一顺时针方向的回流旋涡，随着时间的增加，该回流区的范围逐渐增大；当水库中的水体产生的阻力阻碍其向前进一步运动时，浑水开始下沉（T=1.6时刻），在该回流右侧的自由表面附近的区域内，将出现一逆时针方向的环流；从T=2.4和T=3.2的速度矢量图我们还可看出，随着时间的增加，顺时针回流因受左侧边界的限制，发展逐渐减慢，而右侧环流因异重流沿底部继续地向前运动，异重流与水库水体之间剪切力逐渐发展，范围将逐渐增大，水库水体表面将有倒流现象产生。图6是不同时刻下的相对有效密度等值线图（Γ/Γ0）：当T=0.8时，进口附近区域基本处于初始动量作用的控制，相对有效密度值沿水平方向逐渐延伸，随着时间的增加，进口动量作用越来越小，因密度差而产生的负浮力作用却逐渐增强，Γ/Γ0受负浮力和回流的影响而又逐渐降落，当T=3.2时，Γ/Γ0贴到底部，在底部形成一股底流，直此异重流形成了。 |
|
|
|
|
|
图6 异重流潜入过程中密度场随时间的变化(Fre=1.58) |
图7 T=3.2时刻的水平速度分量沿程分布图(Fre=1.58) |
4.3 潜入前后浑水流速的沿程变化
图7是本文依据数值计算结果，点绘的Fre=1.58的异重流潜入时，在T=3.2时刻的水平速度分量的沿程分布，从该图可以看出，异重流在潜入的过程中，水平速度分布沿程发生了根本的改变：明流最大流速点一般在水面处，潜入后，流速分布型式变化，且最大流速点的相对位置下移，这与文献[8]的分析相同。
5 结论
本文采用混合有限分析法求解了二维Navier Stokes方程和浓度方程，并根据数值计算的结果，对异重流的潜入运动规律进行了分析，结果表明：(1)本文所建立的数学模型能够反映异重流潜入流动的流动特性和内在机制，(2)通过对数值计算结果的讨论和分析可知：①进口密度弗汝德数Fre对异重流潜入形态的影响很大，当Fre&<0.6时，为卷流型，当Fre&>1.5时，为射流型，当0.6≤Fre≤1.5时，流动为浮力射流型；②异重流在潜入的过程中，沿程垂线流速分布发生了根本的改变：最大流速点的相对位置下移。
致谢：本文得到詹义正教授、赵云副教授的指导。
参 考 文 献
[1] 中国学会泥沙专业委员会。泥沙手册。第八章。中国环境科学出版社，1992年。
[2] 钱宁，万兆惠。泥沙运动力学。科学出版社，1983年。第457～460页。
[3] T. H. Yoon， Density Current by Negative Buoyant Inflow in Stagnant， International Association for Hydraulic research， Kyoto Japan， July 20-22，1988， pp.73～80.
[4] 方春明，韩其为，何明民。异重流潜入条件分析及立面二维数值模拟。泥沙研究， 1997，(4).第68～74页。
[5] Evans G and Paolullii S.， The Thermoconvective Instability of Poiseuille Flow Heated from below a Proposed Benchmark Solution For Open Boundary Flows， Int. J. Numer. Methods in Fluids， pp.1001-1013， 1990.
[6] 陶文铨。数值传热学。西安交通大学出版社，1988年。
[7] 李炜。粘性流动的混合有限分析法。研究生教材，武汉电力大学，1994年。
[8] 姚鹏，王兴奎。异重流潜入规律研究。学报，1996，(8).第77～83页。
转贴于