二级沙波水流特性的大涡模拟

董伟良，诸裕良，马殿光，徐俊锋

（1.河海大学港口海岸与近海工程学院，南京210098；2.交通运输部天津水运工程科学研究所，天津300456）

二级沙波水流特性的大涡模拟

董伟良1，诸裕良1，马殿光2，徐俊锋2

（1.河海大学港口海岸与近海工程学院，南京210098；2.交通运输部天津水运工程科学研究所，天津300456）

采用FLuent软件平台大涡模拟方法（LES）对二级沙波水流特性进行数值模拟，进出口边界采用周期性条件以模拟一系列沙波。通过与水槽试验数据对比发现，模拟计算结果与试验数据吻合较好，在此基础上对3种不同水平位置处的二级沙波进行模拟。计算结果表明：当二级沙波距一级沙波波峰较远时，二级沙波对一级沙波起遮掩作用，使得一级沙波回流区范围减小，垂向负流速区更加接近一级沙波波峰，同时二级沙波背流面的水流会对一级沙波迎流面形成侵蚀，进而减小一级沙波尺度；当一、二级沙波接近融合或融合时，两者之间的非线性相互作用使得回流区范围扩大，紊动增强。Q准则等值面显示，二次沙波的叠加增强了沙波涡旋结构的尺度和强度。

二级沙波；大涡模拟；水平位置；Q准则

在河流和海洋环境中，因底床形态不断发展演变，沙波地形也处于不断叠加、融合的过程中，有效地保存了河床底部沉积物的形态［1］。Fernandez［2］试验结果表明：二级沙波沿一级沙波迎流面从波谷移向波峰时，波谷产生的非线性作用使得沙波水流流态处于不断的演变之中。数值模拟方法是研究沙波水流条件及泥沙运动规律的主要手段之一，目前常采用k⁃ε紊流模型，如：Mendoza［3］、Johns［4］、Yoon［5］等人成功模拟了不同条件下沙波水流特性，但模拟结果与实测数据在紊动较强的波谷差距较大；为了克服k⁃ε紊流模型只能模拟时均值及偏差较大的缺点，Yue［6］、Stoesser［7］、Grigoriadis［8］、Xie［9］等选择利用大涡数值模型（LES）来进行研究，相对而言，LES模型可以更为精细地刻画沙波水流条件。

然而，目前沙波水流数值模拟研究主要针对一级沙波形态，较少对二级沙波影响下的沙波水流特性进行研究。Frias等［10］曾对Fernandez［2］的试验过程进行模拟，主要对沙波接近融合和融合2种状态进行研究，而对于二级沙波远离一级沙波波峰的情况暂未涉及。本文拟采用LES模型对二级沙波上的水流特性进行模拟，通过与已有的物理试验资料进行对比，分析水平位置对沙波水流流态的影响，以此研究沙波叠加、融合的演变机制。

1 数学模型

1.1基本方程

大涡模拟的思想是在流动区域内对N⁃S方程进行网格过滤，从而得到较大尺度漩涡的基本方程组。

1.2边界条件

当水流流经一系列沙波时，各个沙波上水流结构存在周期性变化规律，单个沙波周期进出口水流特性基本相同。在保证计算准确性的基础上，为了减少模型网格数和计算时间，从而取为周期性边界条件，以此可以保证在单位周期沙波长度的模拟情况下，沙波水流充分发展为湍流，其中周期性条件为

式中：Φ为纵向流速u、横向流速v、垂线流速w和稳定动能k等水力参数；λ为沙波波长。

由于近壁湍流的脉动能量主要集中在小尺度涡，从而给粗网格的大涡模拟带来困难，为使Smagorinsky⁃Lilly亚网格模型很好地适应近壁面的湍流边界层，采用壁面函数法处理近壁湍流，即用半经验公式将自由流中的紊流与壁面附件的流动链接起来。在固壁上应用无滑移条件，计算域的上边界采用对称边界条件：∂/∂z=0和w=0。

1.3网格生成

计算域由二级沙波叠加在一级沙波迎流面上组成。使用FLUENT软件前处理程序Gambit生成计算区域几何体，再进行网格划分，得到如图1所示的四面体网格单元。模型计算域长为1.6 m，宽为1.0 m，高为0.372 m，采用结构网格进行划分，其中沿x、y和z轴3个方向分别为160、50和70个网格（如图1所示，其中网格每隔3个显示1个），相邻网格比不超过1.05。

图1 模型计算范围及网格剖分Fig.1Calculation scope and computational grid

1.4计算模型

本研究共建4个具有不同沙波形态的数值模型，其概念模型如图2所示。一级沙波采用Mielo和Ruiter［11］试验时的沙波尺度；二级沙波叠加在一级沙波迎流面上，波长0.5 m，波高0.05 m，这与Fernandez［2］试验组合相似。试验过程中，二级沙波处于图2所示的3种不同水平位置（模型b、c和d），二级沙波波谷距一级沙波波峰的长度lS分别为0.5 m、0.16 m和0 m，其中模型a、c和d与Fernandez［2］和Frias［10］试验组合一致，模型b主要针对二级沙波离一级沙波波峰较远时的水流特性进行研究。试验过程中，平均流速uˉ为0.52 m/s，雷诺数Re为171 000，弗汝德数Fr为0.29。

图2 4种床面形态示意图Fig.2Four morphology categories of bed surface

2 计算结果与分析

2.1模型验证

为了验证模型的准确性，选取Mielo和Ruiter［11］的试验数据进行验证。沙波波长1.6 m，波高0.08 m，如图1所示，沙波表面泥沙中值粒径D50=1.6 mm。本文选取T6组次试验数据来对模型进行验证，其中流量Q=0.257 m3/s，水深d=0.334 m，弗汝德数Fr=0.29。

图3显示了不同位置处纵向流速垂线分布情况，其中各测点距前一个波峰分别为0.13 m、0.48 m、0.82 m和1.27 m。水流从波峰流向波谷，再流向波峰，纵向流速剖面因背流面回流区影响垂线分布不均匀，后因迎流面坡度作用水流加速而产生次生流，使得纵向流速垂线分布又趋于均匀［12］。总体而言，LES计算结果和实测值吻合较好，只是在位置0.48 m和0.82 m处，底部流速稍大于实测值。

图3 纵向平均流速垂线分布计算结果与实测值对比Fig.3Comparison of calculation results with measured values of mean longitudinal velocity distribution along vertical

2.2纵向流速

图4为纵向平均流速等值线分布图，流速随过流断面减小而增大，最大流速值出现在二级沙波波峰处。当二级沙波从波谷移向波峰，u=0.6 m/s等值面不断扩大，甚至出现u=0.7 m/s等值面。值得注意的是：当水流流经沙波时，因沙波地形影响，波谷纵向流速u＜0 m/s。以u=0 m/s作为u＜0 m/s的界限，由此可以观察到，当二级沙波距一级沙波波峰较远时（模型b），纵向平均流速u＜0 m/s范围小于模型a；而当一、二级沙波接近融合和融合时（模型c、d），随着二级沙波越接近一级波峰u＜0 m/s范围越大。分析认为：当二级沙波达到一定尺寸后，在其背流面会形成较为明显的回流区，此时因二级沙波距一级波峰较远，一、二级沙波回流没有产生相互作用，二级沙波对一级沙波起到遮掩作用，减小了一级沙波回流区的强度和范围。模型b与Schwammle等［13］提出新月形沙波孤立波现象时的水流条件相似，当二级沙波尺度较大时，其背流面的水流会对一级沙波形成侵蚀，同时二级沙波阻碍了一级沙波的来沙，致使一级沙波尺寸不断减少，二级沙波尺寸不断增加，从而形成孤立波现象。

2.3垂向流速

图4 纵向平均流速等值线分布图Fig.4Distribution of mean longitudinal velocity isoline

图5给出了4种模型的垂向流速等值线分布图，由图5可见，尽管二级沙波的叠加改变了一级沙波迎流面上的分布情况，但一般分布规律并没有改变：因沙波的迎、背流面地形影响，在迎流面会出现垂向流速w＞0 m/s区域（标记为U），背流面垂向流速w＜0 m/s区域（标记为D）。垂向负流速区D因u＞0 m/s、w＜0 m/s，水流冲击床面，泥沙运动较为剧烈，床面侵蚀也剖为严重。

图5 垂向平均流速等值线分布图Fig.5Distribution of mean vertical velocity isoline

因二级沙波位置影响，一级沙波背流面垂向负流速区D发生明显改变。模型b中，一、二级沙波的垂向负流速区D没有融合，其中以w=0 m/s和0.04 m/s等值线作为参考，从图5中可以看出一级沙波背流面后的垂向负流速区D相对于模型a明显减小。研究认为，垂向负流速区D的减小主要由以下2个因素引起：（1）二级沙波对一级沙波起遮掩作用；（2）二级沙波叠加在一级沙波迎流面上时，二级沙波上形成的垂向正流速区U在一定程度上限制了垂向负流速区D的发展。与Fernandez［2］试验的试验结果相同，在二级沙波愈加靠近一级沙波波峰的情形下，一、二级沙波背流面垂向负流速区D相融合，同时由图5可以看出，D范围也愈大，垂向流速w=-0.06 m/s等值面也离一级沙波波峰愈远，甚至还出现了w=-0.08 m/s的垂向流速。

因一级沙波波陡λ=1/20小于二级沙波波陡λ=1/10，将垂向流速w=0.02 m/s作为垂向正流速区U的界限，可以看出：模型a的垂向正流速区U相对比较扁平，模型b因空间有限，受上游垂向负流速区D的挤压，垂向正流速区U最小；随着二级沙波愈加靠近一级沙波波峰，和垂向负流速区D一样，垂向正流速区U范围也逐渐增加。

2.4涡旋结构

剪切层Kelvin⁃Helmholtz不稳定性是沙波产生涡旋的主要原因［7、9］，图6为某一时刻4种不同床面形态下横向位置y=0.5 m截面上横向涡量ωy的分布图，其中在一、二级沙波背流面均形成了较明显的大尺度涡旋。从图6中可以看出，模型a中横向涡量ωy主要贴附底床，上层水体中ωy值普遍偏小；模型b中，一、二级沙波产生的横向涡量ωy并没有相互作用，在二级沙波背流面可以明显观测到二级沙波产生的横向涡量ωy，较强的涡旋会对床面造成侵蚀；模型c、d中，一、二级沙波产生的横向涡量ωy相互作用，其中模型d中，在上层水体中也可以观察到较大的横向涡量ωy。

图6 某一时刻横向涡量ωy分布图Fig.6Distribution of transverse vorticityωyat a specific time

为了更清楚地显示流场中的涡旋特征及破碎过程，采用Q准则［14］方法进行识别。图7显示了Q=9对应的涡旋结构。从图7中可知，涡旋结构主要集中在波谷附近出现，根据前人研究分类，沙波水体中涡旋结构一般可分为：横向涡、发卡（马蹄）涡和纵向涡3类［8-9］。由于剪切层Kelvin⁃Helmholtz不稳定性，在一、二级沙波波峰处首先会出现横向涡，随着向下游发展，横向涡逐渐演变成发卡涡。随着流场进一步向下游发展，发卡涡逐渐向斜上方演变，由于上层水体流速大，下层流速小，涡旋结构被不断拉长扭曲，其中一部分到达水面形成水面翻滚现象，一部分演变成纵向涡进而破碎成随机分布的小蜗，并逐渐消失。

当只有一级沙波时（模型a），与横向涡量ωy显示结果相同，涡旋结构主要集中在波峰线以下水体。当二级沙波叠加时，背流面涡旋结构尺度变大，其中在二级波峰顶处可以明显观察到涡旋结构B，这主要是由于二级沙波产生的涡旋结构在上部相对较为平静的水体中，随着流场的发展，能够较好地保持结构形态而不破碎。图7显示，当二级沙波离一级沙波波峰较远时（模型b）与接近融合和融合（模型c、d）状态不同，二级沙波产生的涡旋结构没有能够继续发展，部分已经破碎，同时在上层水体中也没出现模型c、d中大尺度的涡旋结构K。

图7 某一时刻涡结构分布图（Q=9）Fig.7Distribution of vortex structure at a specific time（Q=9）

3 结论

本文利用LES模型对不同水平位置处的二级沙波进行数值模拟。比较本文的计算结果和实测数据，可以看出LES可以准确刻画沙波上的水流特性。二级沙波从波谷移向波峰时，其对一级沙波的作用经历了从遮掩侵蚀到非线性相互作用的演变：当二级沙波离一级波峰较远时，二级沙波形态对下游一级沙波起遮掩作用，使得一级沙波回流区面积缩小，垂向负流速区域更加靠近一级波峰，同时二级沙波背流面的水流会对一级沙波迎流面形成侵蚀，甚至会减小一级沙波尺度；当一、二级沙波接近融合或融合时，两者非线性作用使得背流面回流区面积扩大，紊动增强。Q准则等值面显示，二次沙波的叠加使得沙波水体中涡旋结构尺度和强度均匀增加，因上层水体紊动较小，二级沙波产生的涡旋结构更易保持形态而不破碎。

［1］Best J，Blois G，Barros J，et al.The dynamics of bedform amalgamation：new insights from a very thin flume［C］//Fourth Internation⁃al Conference on Marine and River Dune Dynamics.Belgium：Bruges，2013：29-34.

［2］Fernandez R，Best J，López F.Mean flow，turbulence structure，and bed form superimposition across the ripple⁃dune transition［J］. Water Resources Research，2006，42（5）：72-88.

［3］Mendoza C，Wen S H.Investigation of turbulent flow over dunes［J］.Journal of Hydraulic Engineering，1990，116（4）：459-477.

［4］Johns B，Soulsby R L，Xing J X.A comparison of numerical model experiments of free surface flow over topography with flume and field observations［J］.Journal of Hydraulic Research，1993，31（2）：215-228.

［5］Yoon J Y，Patel V C.Numerical Model of Turbulent Flow over Sand Dune［J］.Journal of Hydraulic Engineering，1996，122（1）：10-18.

［6］Yue W，Lin C L，Patel V C.Numerical investigations of turbulent free surface flows using level set method and large eddy simula⁃tion［R］.Iowa City：Iowa Inst.of Hydraul.Res，2003.

［7］Stoesser T，Braun C，Garcia-Villalba M，et al.Turbulence structures in flow over two⁃dimensional dunes［J］.Journal of Hydraulic Engineering，2008，134（1）：42-55.

［8］Grigoriadis D G E，Balaras E，Dimas A A.Large⁃eddy simulations of unidirectional water flow over dunes［J］.Journal of Geophysical Research：Earth Surface，2009，114（F2）：91-100.

［9］Xie Z，Lin B，Falconer R A.Turbulence characteristics in free⁃surface flow over two⁃dimensional dunes［J］.Journal of Hydro⁃envi⁃ronment Research，2014（8）：200-209.

［10］Frias C E，Abad J D.Mean and turbulent flow structure during the amalgamation process in fluvial bed forms［J］.Water Resources Research，2013，49（10）：6 548-6 560.

［11］Mierlo M，De Ruiter J C C.Turbulence measurements above artificial dunes［R］.Delft：Delft Hydraulics Lab，1988.

［12］马殿光，董伟良，徐俊锋.沙波迎流面流速分布公式［J］.水科学进展，2015，26（3）：396-403. MA D G，DONG W L，XUN J F.Velocity distribution of nonuniform flow on the stoss side over dunes［J］.Advances in Water Sci⁃ence，2015，26（3）：396-403.

［13］Schwammle V，Hermann H J.Solitary wave behaviour of sand dunes［J］.Nature，2003，426（11）：619-620.

［14］Jeong J，Hussain F.On the identification of a vortex［J］.Journal of Fluid Mechanics，1995，285（4）：69-94.

Large⁃eddy simulation of flow characteristics over secondary dunes

DONG Wei⁃liang1,ZHU Yu⁃liang1,MA Dian⁃guang2,XU Jun⁃feng2
（1.College of Harbor,Coastal and Offshore Engineering，Hohai University，Nanjing 210098，China;2.Tianjin Research Institute for Water Transport Engineering,M.O.T.,Tianjin 300456,China）

Based on the Fluent software,the large eddy simulation(LES)technique was used to model the flow characteristics over secondary dunes.Periodic boundary condition was imposed at the inlet and outlet to simulate a series of dunes.Good agreement was obtained between the simulation and experimental data,and secondary dunes at three different positions were simulated.The results show that the secondary dune has a shelter effect on the pri⁃mary dune when far from the crest of the primary dune,which reduces the range of recirculation zone and makes ver⁃tical negative velocity zone more close to the primary dune crest.Additionally,the current at the lee side of second⁃ary dunes can cause erosion at the upstream side of the primary dune,thus diminish the dimension of the primary dunes.However,recirculation zone becomes large and turbulence intensities increase due to nonlinear interaction between current and dunes when secondary dunes approach and merge with the primary dune.The size and strength of vortex structures,visualized by the isosurface of Q⁃criterion,were strengthened by the secondary dune.

secondary dune;large eddy simulation;horizontal position;Q⁃criterion

TV 135；O 35

A

1005-8443（2016）02-0154-05

2015-07-29；

2015-11-30

中央级公益性科研院所基本科研业务费专项资金项目（TKS150102，TKS130105）

董伟良(1990-)，男，江苏省人，硕士研究生，主要从事河流水力学、港口航道研究。

Biography：DONG Wei⁃liang（1990-），male，master student.