kinoshita派生弯道水流三维数值模拟

白玉川，徐国强

（1.天津大学 水利工程仿真与安全国家重点实验室，天津 300072；2.天津大学 河流海岸工程泥沙研究所，天津 300072）

kinoshita派生弯道水流三维数值模拟

白玉川1,2，徐国强1,2

（1.天津大学 水利工程仿真与安全国家重点实验室，天津 300072；2.天津大学 河流海岸工程泥沙研究所，天津 300072）

本文采用雷诺方程和雷诺应力紊流模型相结合建立起弯道的三维数学模型，并结合shukry弯道实验结果对模型进行验证。利用此模型对 kinoshita型非对称弯道进行模拟研究。结果表明：非对称型弯道的高流速区比对称型弯道更加辐聚，水流动力轴线更加的不顺滑，主流区在弯顶的顶冲作用更强；对称型弯道的断面环流强度不断衰减，而非对称型弯道则不同；随着偏斜度的增大，断面环流强度在弯道前半段不断减小，后半段则相反。

非对称型弯道；三维数学模型；kinoshita曲线；偏斜度

引 言

天然弯曲河流的平面形态多呈现出非对称的平面形态，这对弯道水流有着重要的影响。现有的研究多是针对对称型弯道进行研究，而对于非对称型弯道的研究较少。kinoshita派生曲线与非对称型河流的平面形态吻合较好，本文根据kinoshita派生曲线建立多组非对称型弯道，并采用雷诺应力模型对其进行数值模拟，研究非对称型弯道的水流运动特性。

1 非对称型弯道的三维数学模型

1.1 控制方程

描述不可压缩流体的控制方程是雷诺方程（Reynolds equation），包括连续性方程和动量方程。为了使雷诺方程封闭，本文采用Reynolds应力紊流模型。

连续性方程：

动量方程：

Reynolds应力运输方程：

k方程：

ε方程：

12C1ε=1.44，C2ε=1.92，Cμ=0.09，σk=0.82，σε=1.0。

1.2 模拟自由表面的VOF方法

VOF方法的基本思想：用体积率函数表示不同流体在计算单元中所占的比例。对于两相流来说，分别定义aw(x,y,z,t)和aa( x,y,z,t)表示水和空气在计算单元中的体积比。在任何一个计算单元中水和空气的体积比之和均为 1。因此任何一个单元会有全是水、充满空气或水气都有三种情况。

通过求解下面的连续性方程来跟踪水气界面：

引入VOF后，自由表面处的密度ρ和动力粘度μ不在是常数了，而是体积分数的函数了。它们可以由下式计算得到：

式中：μw和ρw分别表示水的动力粘度和密度，μa和ρa分别表示空气的动力粘度和密度。

1.3 边界条件与数值计算方法

本文采用有限体积法对方程进行离散，为了提高计算精度，网格采用结构化六面体网格。网格首先在s-n方向剖分，网格垂直于弯道中心线，然后在垂向（z方向）上加密生成结构化六面体网格，网格三维效果如图1。对于压力和速度场采用PISO算法进行耦合计算。

图1 弯道网格剖分三维效果

弯道进口分为水速度进口和空气压力进口，水速度进口采用常数分布，弯道出口采用压力出口，弯道底部和侧壁均采用固定边壁，粘性底层采用壁面函数处理。

1.4 模型的验证

采用Shukry[1]的弯道实验对模型进行验证。图2为模型弯道纵向垂线平均速度的实测值（图2a）和计算值（图2b）的对比图。从图2中可以发现，实测值和计算值的流速分布吻合较好，高流速区在弯顶凸岸附近，弯顶下游低流速区也吻合较好。本文采用的模型适合于弯道水流的模拟，本文采用上述模型对非对称型弯道进行模拟，研究非对称型弯道的水流特性。

图2 等值线实测值与计算值对比

2 非对称弯道的数值模拟

kinoshita派生曲线[2]考虑了河流的偏斜度和丰盈度，描述弯道的平面参数更加的完备，可以通过调整最大偏角、弯道的波长、丰盈度系数、偏斜度系数四个参数来描述弯曲河流的平面形态，kinoshita派生曲线能够很好的吻合天然弯曲河流的平面形态。

式中：Js为偏斜度系数；Jf为丰盈度系数；L为弯道的波长；s为弯道水流方向的坐标。

偏斜度系数对弯道的非对称特性有着重要的影响，因此本文通过改变偏斜度系数，建立不同的非对称型弯道，探究其水流特性。具体的模拟工况如表1，图3为各个工况下弯道中心线，图中可以看出，S1为对称型弯道，S2、S3为非对称型弯道。

表1 不同工况的参数设置

图3 各个工况下弯道中心线

2.1 纵向流速

纵向流速对弯道水流有着重要作用，为了分析非对称弯道下的纵向流速，分别绘制各个工况弯顶断面的纵向流速。

图4 各个工况下弯顶断面纵向速度

从图中可以看出，纵向速度在横向上存在明显的坡降，高流速区（2倍的平均速度）主要集中在靠近凸岸的0.1倍河宽范围内，凹岸的速度明显小于凸岸速度。S1工况下弯顶处断面的高流速中心主要集中 0.4～0.6倍的水深处，并且比较对称。S2、S3弯顶断面的高流速中心则是不断的下移，靠近弯道水槽底部。非对称型弯道的弯顶高流速区中心比对称型弯道更加靠近水槽底部。可以推测，随着偏斜度的增大，水流会不断的侵蚀凸岸的底部，可能会在此形成崩岸，促使河湾由曲变直，使河湾加快消亡。

2.2 水流动力轴线

水流的纵向垂线平均速度是弯道水流一个重要特征，可以整体反映水流的特性。

图5为各个工况下的纵向垂线平均速度等值线图，图中实线为水流的动力轴线。从图中可以看出对称型弯道和非对称型弯道均会在弯顶凸岸附近形成高流速区，而在凹岸形成低流速区。随着偏斜度的增大，高流速区在弯顶凸岸越来越辐聚，流速也在不断的增大，而紧靠弯道凸岸下游流速会较小，当偏斜度达到一定程度时，此区域会形成回流。

S1工况下的水流动力轴线相对比较顺滑，形态上比较对称。随着偏斜度的增大，水流动力轴线越来越不顺滑，与弯顶凸岸的夹角越来越大，此时弯道主流与弯顶的顶冲作用较强。可以推测，回流区的形成会使泥沙在凸岸进一步的淤积，有利于凸岸边滩的发展，进一步的促进弯道的发展，当弯道发展到一定的程度时，高流速区域直接切割凸岸边滩，形成裁弯。

图5 纵向垂线平均速度及水流动力轴线

2.3 断面环流强度

将弯道沿弯道中心线均分成31个断面，S1工况断面布置如图6。横向环流强度没有明确的定义，由于弯道中的水流运动是一种螺旋流，引入气象学中的螺旋度概念表示环流强度。螺旋度（H）大小反映的是旋转与沿旋转轴方向的强弱程度[3]。螺旋度的计算公式如下：

为了更好的分析环流强度沿程的变化，将环流强度的绝对值在断面上进行积分，得到断面的环流强度（DH）。

图6 断面布置

图7 断面环流强度沿程分布

从图中可以看出，S1工况下，断面环流强度会出现三个峰值且峰值是不断衰减，沿程呈现出一种类阻尼现象，而非对称型弯道则不同。对比S1、S2和 S3可以看出，随着偏斜度的增加，弯道断面环流强度在弯道的前半段是不断地减小，而在弯顶下游，尤其是出口处，断面环流强度逐渐地变大。

2.4 床面切应力

弯道的床面切应力对于泥沙的输运有着重要的作用，研究弯道的床面切应力分布可以预测泥沙在弯道中的冲淤情况。床面切应力按照下面的公式计算[4]：

式中：Uw为流速；y为到床面的距离；ks为壁面粗糙度厚度。

图8 床面切应力分布

从图中可以看出，S1弯道的切应力的范围为0.1至0.85之间，在弯顶处切应力在横向方向上存在明显的坡降，床面切应力在弯顶凸岸较大，而在弯顶的凹岸切应力则相对较小，切应力在弯道的两侧基本呈现出对称型分布。S3弯道在进口右岸的切应力较低，切应力在弯道的上半段相对于明显小于下半段，在弯顶凸岸的下游存在着较小的切应力区域。非对称型弯道，切应力在弯顶凸岸处明显比S1弯顶大。在弯道的前半段，S3切应力在横向上存在明显的坡降，而 S1弯道虽然存在着坡降，但没有S3的坡降大。对比S1、S2、S3可以看出，随着弯道偏斜度的增大，弯顶处凸岸的切应力不断的增大，并且低切应力的区域范围也在不断增大。对比图5和图8，可以看出，弯道床面切应力分布与纵向垂线平均流速的分布基本一致，弯顶的高流速区域切应力相对也较大，而低流速区域切应力相对也较小。

3 结 论

本文采用kinoshita派生曲线建立多组非对称型弯道，并采用数值模拟的方法对非对称型弯道进行研究，根据模拟结果分析，得出以下结论：

1）对称型弯道和非对称型弯道在弯顶凸岸处均存在高流速区，随着偏斜度的增加，高流速区在弯顶凸岸会更加辐聚，非对称型弯道的高流速区比对称型弯道更加集中。此外水流动力轴线越来越不顺滑，主流区与弯顶的顶冲作用较强。可以推测，随着偏斜度的增大，弯顶凸岸附近的高流速区，可能切割凸岸边滩，形成切滩型裁弯。

2）随着偏斜度的增加，弯道断面环流强度在弯道的前半段是不断地减小，而在弯道的后半段，断面环流强度逐渐地变大。

3）弯道床面切应力分布与纵向垂线平均流速的分布基本一致。对称型弯道的床面切应力分布基本对称，而非对称弯道的床面切应力则不是，随着偏斜度的增大，弯道前半段的切应力小于弯道后半段的切应力。

[1]Shukry A.Flow around bends in an open flume[J].Transactions of the American Society of Civil Engineers,1950,115:751-779．

[2]Kinoshita R．Investigation of channel deformation in Ishikari River[R]．Report to the Bureau of Resources,1961:1-174．

[3]陆慧娟,高守亭.螺旋度及螺旋度方程的讨论[J].气象学报,2003,61(6):684-691.

[4]Olsen N R B.A three-dimensional numerical model for simulation of sediment movements in water intakes with multiblock option[J].User’s Manual,The Norwegian University of Science and Technology,Trondheim,2006.

3D Numerical Simulation of Water Current along Kinoshita-generated Bend

Bai Yuchuan1,2，Xu Guoqiang1,2
(1.State Key Laboratory of Hydraulic Engineering Simulation and Safety,Tianjin University,Tianjin 300072,China；2.Institude of Sediment on River and Coast Engineering,Tianjin University,Tianjin 300072,China)

Reynolds equation and Reynolds stress turbulence model are combined to build a 3D mathematical model of bend,and the model is verified according to shukry bend testing results.Then this model is used to study kinoshita asymmetric bend.The results show that high flow velocity zone of asymmetric bend is more concentrated than that of symmetric bend,water current dynamic axis of the former one is not smooth as the latter one,and the impact force of main flow is greater at the apex of bend.The intensity of circulating flow on section of symmetric bend weakens continuously.With the increase of skewness coefficient,the intensity of circulating flow on the section of asymmetric bend decreases in the first half of bend and increases in the second half of bend..

asymmetric bend; 3D mathematical model; Kinoshita-generated curve; skewness

U611

：A

：1004-9592（2016）06-0007-05

10.16403/j.cnki.ggjs20160602

2015-12-13

自然科学基金资助项目（51279124）

白玉川（1967-），男，教授，主要从事泥沙运动力学及海岸动力学研究。