细颗粒泥沙净冲刷和输移的大涡模拟研究

白 静方红卫何国建谢崇宝高 虹

∗(中国灌溉排水发展中心，北京100054)

†(清华大学水利水电工程系，水沙科学与水利水电工程国家重点实验室，北京100084)

∗∗(北京中灌绿源国际咨询有限公司，北京100054)

细颗粒泥沙净冲刷和输移的大涡模拟研究

白 静∗方红卫†,1)何国建†谢崇宝∗高 虹∗∗

∗(中国灌溉排水发展中心，北京100054)

†(清华大学水利水电工程系，水沙科学与水利水电工程国家重点实验室，北京100084)

∗∗(北京中灌绿源国际咨询有限公司，北京100054)

在传统水沙输移数值模拟研究中一般采用雷诺时均模拟技术(Reynolds-averaged simulation,RANS).与RANS相比，大涡模拟技术(large eddy simulation,LES)能够更加精确反映细部流动结构，计算机的发展使得采用LES探讨水流和泥沙运动规律成为可能.本文尝试给出净冲刷条件下悬沙计算的边界条件，采用动态亚格子模式对循环槽道和长槽道中的水流运动和泥沙输移进行了三维大涡模拟研究.利用直接数值模拟(direct numerical simulation,DNS)结果对LES模型进行了率定，计算结果符合良好，在此基础上初步探讨了泥沙浓度、湍动强度和湍动通量等的分布特征.结果表明，净冲刷条件下输沙平衡时泥沙浓度符合Rouse公式分布，单向流动中泥沙浓度沿着流向逐渐增大.泥沙浓度湍动强度和湍动通量都在近底部达到最大值，沿着垂向迅速减小.湍动黏性系数和扩散系数基本上在水深中间处达到最大.湍动Schmidt数沿着水深方向不是常数，在近底部和自由水面附近较大，在水深中间处较小.

细颗粒泥沙,净冲刷,大涡模拟,泥沙输移

引言

泥沙输移一直是水力学及河流动力学领域研究的热门问题.淤积和冲刷是不平衡输沙的重要形式.不平衡输沙是河床变形和水库淤积的根源.细颗粒泥沙的输移一直受到广泛的关注[1-2].细颗粒泥沙在水体污染、水体富营养化和水生生态系统健康等方面有重要的影响.来自农业、工业和生活等污水进入河道后，以细颗粒泥沙为载体的生物膜生长过程，改变了泥沙的形态特征、沉降和起动特性[3].细颗粒泥沙的吸附解析作用和絮凝过程，对河道和水库中污染物的迁移转化[4]和重金属的迁移[5]有复杂的影响.泥沙输移的计算方法主要有两种：欧拉方法[6-7]和拉格朗日方法，欧拉方法的计算量大，只能模拟有限泥沙颗粒的运动.在泥沙输移和河床变形模拟中，泥沙颗粒的数量巨大，因此应当采用拉格朗日方法进行模拟.目前绝大多数的泥沙运动和水流模拟成果主要采用雷诺时均模拟(Reynolds-averaged simulation,RANS)完成[8-10]，Keylock等[11]认为随着计算机技术水平的提高，使采用大涡模拟(large eddy simulation,LES)追踪和研究湍流相干结构的变化成为可能，由于湍流相干结构对泥沙起动和湍流输运有极其重要的作用，LES在泥沙输移和河床演变研究中的潜力和优势也逐渐表现出来.根据前人的研究，对于同样的算例,LES的计算网格量为三维雷诺时均模拟的数十倍乃至数百倍[12]，计算时间和计算量随着计算区域长度的增大而增加.

国内外的学者针对LES中泥沙计算边界、泥沙起动的机理、浓度分布特征和沙纹的形成过程进行了探讨.在定床计算中，泥沙颗粒一般较细，一般只涉及到悬移质.Zedler等[13-14]假定底部泥沙浓度的梯度等于pick-up函数与扩散系数 ε的商，即∂s/∂z=-Pk/ε，先后模拟了单向流和振荡流中泥沙起动过程，在单向流中得到了泥沙浓度随时间的变化过程，得到了振荡流中流速、压强以及泥沙浓度分布.其研究认为泥沙的起动原因有两种：一种是泥沙在涡体形成过程中被带起，一种由于床面切应力较高，泥沙直接起动.Chou等[15]利用动态混合模型和LES计算了泥沙在明渠中起动直至平衡的全过程，平衡时的泥沙浓度与Rouse公式符合较好，认为通过床面进入水体的泥沙通量主要由流向涡决定.Widera等[16-17]利用LES计算平衡时泥沙在平底和有沙波的条件下槽道水流中的泥沙分布，主要探讨了湍动Schimidt数的取值问题，指出了在RANS中将不同水深处湍动Schmidt数取为常数的做法是不可靠的；在有沙波存在的明渠中，底部泥沙浓度受地形的影响较大，在远离床底处，泥沙浓度与Rouse公式符合很好.靠近床面处，湍动Schmidt数对地形不是很敏感.

在前面所提及的利用LES对泥沙分布与输移的研究中，Zelder等[13-14]采用的泥沙边界条件，在扩散系数ε较小时，容易造成泥沙浓度s在底部的梯度过大.Chou等[15]引入泥沙沉降通量Db和起动通量Eb的概念，Eb通过pick-up函数计算得到，近底处的泥沙浓度sb通过竖直方向上邻近的三点浓度插值得到，这种泥沙边界的处理方法，在计算时没有用到参考高度的概念，比较简单，但是Eb实际上是全沙通量，并不准确.由于流向采用周期性边界条件，只能给出泥沙浓度和床面变形随时间的变化过程，不能给出空间变化过程.Widera等[16-17]在水平方向上采用周期性边界条件，在底部泥沙计算中采用零通量的边界条件，使得计算结果只能反映泥沙输沙平衡后的最终状态，不能反映泥沙随着时间和空间的变化情况.

泥沙可以分为悬移质和推移质，推移质集中在床面附近.与泥沙的净淤积边界条件相比，净冲刷条件相对比较复杂[18].为了解决净冲刷过程中悬沙的起算高度问题，在计算中保证泥沙和水流计算网格的一致性，本文中将借鉴RANS的处理方式，在悬沙计算中引入参考高度概念，建立细颗粒泥沙净冲刷条件下的近底边界条件，尝试给出输沙平衡条件下和单向流动中的水流运动和泥沙输移规律，对其中水流和泥沙的时均特征、湍动特性进行分析和讨论，进一步统计湍动黏性系数和扩散系数的分布，给出湍动Schmidt数的分布规律.

1 控制方程

水流计算的控制方程为

其中“¯”表示过滤后的变量，ui(i=1,2,3)表示i方向上的无量纲速度；υ是雷诺数Re的倒数；p是无量纲压强；对控制方程中的非线性对流项过滤产生的亚格子应力为它代表亚格子尺度运动对可解尺度运动的影响.

其中υSGS为亚格子黏性系数为 Kronecker算子，本文采用Germano等[19]提出的动态模式进行求解亚格子黏性系数υSGS.

悬沙运动的控制方程为

其中s为无量纲的悬沙浓度，ε为悬沙的分子扩散系数，ws为悬沙的无量纲沉速，ε=υ/σ，σ是Schmidt数，本文中取代表了亚格子尺度的水流运动对悬沙的输运作用，ϑi为LES中悬沙控制方程的不封闭项.

其中εSGS为泥沙亚格子尺度的扩散系数，εSGS与亚格子尺度黏性系数有关：εSGS=υSGS/σSGS，σSGS是亚格子尺度的Schmidt数，本文中取1.

2 数值方法与悬沙参考高度

对德国卡尔斯鲁厄大学流体力学研究所开发LESOCC2程序[20-23]进行了修改，加入了悬沙计算模块.在程序中采用有限体积法对方程进行离散，变量布置在非交错曲线网格上，分别采用二阶精度的中央差分和HLPA差分[24]近似水流和泥沙计算控制方程中的对流通量，统一采用中央差分近似扩散通量，时间推进利用二阶精度的3步Runge--Kutta法.在Runge--Kutta法的第3步，采用Stone[25]提出的强隐式迭代(strongly implicit procedure,SIP)方法求解压力修正的泊松方程.

推移质集中在床面表面附近区域，悬移质计算的网格起点也就是参考高度a，应该不低于推移质层的高度.不同的研究学者对参考高度有着不同的研究，Celik等[26]认为：对于没有沙波的槽道，参考高度为a=2/3ks=2d50，其中ks为床面的当量粗糙高度；对于有沙波存在的槽道，参考高度为a=2/3Δ，其中Δ为沙波的高度.van Rijn[27]认为对于没有沙波的槽道，参考高度为a=ks，对于有沙波存在的槽道，参考高度为a=Δ，并在他的研究中a≥0.01H，H为水深.在本文中采用van Rijn提出的方法确定参考高度.在LES中，由于在近床面处的网格较小，因此悬沙浓度从离开床面的第n个网格开始计算，如图1中加粗的网格所示.

图1 计算网格示意图Fig.1 The sketch map of the computational grids

3 计算算例和边界条件

在本文中设置 2个算例，雷诺数Re均为20000，建立相应于流向、展向和垂向的x×y×z的坐标系.算例1为循环算例，计算区域大小为6H×3H×H(H表示水深)，x，y和z向的网格数分别为210，210和126.x和y向采用均匀网格，采用摩阻流速无量纲化后网格尺寸分别为Δx+=28.08和Δy+=14.04，z向采用不均匀网格，从逐渐扩大到在计算中首先计算120个无量纲时间(H/U)，待计算稳定后，统计了300个无量纲时间，在大型并行计算机上采用36个核计算完成.

算例1中，在流向上和展向上采用循环边界条件.水流计算中在自由水面处引入刚盖假定，施加滑移边界条件.泥沙离开或进入自由水面的净通量为零，在自由水面处采用零通量边界条件 (式 (6)).在床面处水流计算中施加无滑移条件，采用Werner等[28]提出的壁函数计算最靠近壁面的第一层网格的流速值

假定LES中近底部瞬时泥沙浓度符合与RANS相同的计算公式，同时为了使得计算不发散，引入εl有

其中sb和sn分别为为近网格底部和中心的泥沙浓度，sb∗为近网格底部的挟沙力，zb和zn分别为近网格底部和中心的高程.εl为大尺度运动扩散系数，表示大涡模拟能够识别的包含在大尺度运动中的扩散系数

其中，-u′w′为zb处的切应力，∂u/∂z为相应位置瞬时流向流速u的z向梯度.σl为大尺度运动Schmidt数，取为1.为了保证计算稳定性，大尺度运动扩散系数εl采用统计平均值

其中〈〉表示统计平均.

算例2为单向流动算例，计算区域大小为48H× 3H×H，分为先导计算域和主计算域，其中先导计算域的大小为6H×3H×H，在先导计算域中只计算水流，不计算泥沙.在算例2中，x,y和z向的网格数分别为1680,210和126.网格尺寸同算例1.在计算中首先计算288个无量纲时间，待计算稳定后，统计了1800个无量纲时间，在大型并行计算机上采用288个核计算完成.

算例2中，水流计算边界条件的设置与算例1相同.在主计算域中计算水流和泥沙，展向和垂向边界条件与算例1相同.流向上的边界条件设置为：在主计算域的进口处，设置进口边界条件，使得在先导计算域中充分发展的湍流顺利进入主计算域；在出口处，泥沙和水流都采用对流边界条件，并且实时校正水流流量，保证流量守恒.详见图2.

图2 算例2计算区域平面示意图Fig.2 The sketch map of the computational domain in case 2

4 模型的率定与验证

利用算例1的水流计算结果对LES模型进行率定和验证.采用算例 1的结果和直接数值模拟(direct numerical simulation,DNS)[29]计算结果进行比对，流向流速和湍动强度与DNS的比较分别见图3(a)和图3(b).在图3(a)中纵坐标为u+，即对流向的平均流速采用摩阻流速u∗进行无量纲化处理有u+=〈u〉/u∗，〈u〉表示对流向流速取统计平均，摩阻流速u∗=0.04915U，U为断面平均流速.横坐标为z+，采用摩阻流速u∗和υ进行无量纲化.在计算中，基于摩阻流速u∗的雷诺数为Re∗=u∗H/υ=983，DNS中Re∗=950.从图中可以看出，LES计算得到的结果与DNS结果符合良好，仅在自由水面处有细微的差别，LES计算得到的u+比DNS结果偏小0.8.不难发现，当z+介于[30，200]时，LES计算得到的u+服从对数律分布，这与以往的研究结果是一致的.

图3 算例1中LES和DNS流向速度和湍动强度比较Fig.3 Comparisons of the stream wise velocity and turbulence intensities between LES and DNS in case 1

图3 算例1中LES和DNS流向速度和湍动强度比较(续)Fig.3 Comparisons of the stream wise velocity and turbulence intensities between LES and DNS in case 1(continued)

在图 3(b)中比较了流向、展向和垂向上的无量纲湍动强度的分布.其中从图中可以发现，在整体上LES计算的湍动强度与DNS结果符合较好，在床面附近，LES计算得到的的最大值为2.96，比DNS偏大0.14.和w′+的最大值分别为1.33和1.04，比DNS分别偏小0.12和0.04.在z+介于[200,800]时，LES计算得到的略微偏小，而与DNS结果基本一致.由于LES中在自由水面处采用刚盖假定，并施加滑移(对称)边界条件，因此水面处的而略微上抬.

图4 无量纲切应力计算结果对比图Fig.4 The comparison of the dimensionless shear-stresses between LES and DNS

5 结果分析

5.1 时均的水流结果分析

由于循环槽道的流速分布等已有了详细的讨论，这里集中分析长槽道算例(算例2)中的水流结果.由于在长槽道的展向y上采用循环边界条件，水流在展向上是均匀的，因此分析采用展向平均的结果.

图5(a)中给出的是长槽道中x=0,10H,20H, 30H,40H处的瞬时流向流速u/U和时均流向流速〈u〉/U沿着水深的分布情况，均做展向平均处理，其中细实线表示瞬时流向流速u/U，粗实线表示时均流向流速〈u〉/U.从图中可以看出，不同位置处时均流向流速〈u〉/U相同，瞬时流向流速u/U沿x轴也差别不大.时均流向流速〈u〉/U沿x方向均匀分布，在水面处取得最大值1.15，采用摩阻流速u∗无量纲化后为23.3，与图3(a)的结果一致.

图5(b)中给出的是长槽道中x=0,10H,20H, 30H,40H处展向平均的湍动强度分布情况，其中从图中可以看出沿着流向各个位置的湍动强度基本重合，湍动强度在床面处为0，随着水深急剧增加，在近底部取得最大值3.0，1.32和1.04，之后随着水深的增加不断减小.在水面处，垂向湍动强度为0，相应地，流向和展向湍动强度略有抬升，在总体上，长槽道中湍动强度的分布与算例1中的湍动强度分布基本相同(图3(b)).在算例2中，长槽道单向流为均匀流动.

图5 算例2中不同位置上流向速度和湍动强度分布Fig.5 The profile of the stream-wise velocities and turbulence intensities at dif f erent locations in case 2

图6中给出的是长槽道中x=0,10H,20H,30H, 40H处的无量纲切应力分布，从图中可以看出，沿程的剪切应力基本一致.在底面和水面处均为0，离开床面处，迅速减小，-0.89为其最小值，之后沿着水深基本上呈现线性增加的趋势，直到水面处归零，与算例1中结果基本相同.

5.2 时均泥沙浓度计算结果分析

图7(a)中给出的是循环槽道中 (算例 1)无量纲的时均泥沙浓度沿着水深的分布，其中横坐标为〈s〉/sb，sb为输沙平衡时近底部的参考泥沙浓度.算例1可以给出输沙平衡时泥沙浓度沿着水深方向的分布情况.从图中可以看出，在输沙平衡时LES计算得到的泥沙浓度符合Rouse公式的分布.

图6 算例2中无量纲切应力的分布Fig.6 The profile of the dimensionless shear stresses at dif f erent locations in case 2

图7 时均泥沙浓度分布Fig.7 The profile of time-averaged sediment concentration

在图7(b)中，从左到右依次为：算例2长槽道中x=0,10H,20H,30H,40H处的展向平均的〈s〉/sb和算例1中输沙平衡时统计得到的〈s〉/sb.从图中可以看出在长槽道中，泥沙浓度始终小于输沙平衡时的泥沙浓度〈s〉/sb，可以判断出长槽道处于净冲刷状态.沿着x方向，随着冲刷的进行，泥沙浓度由于受到水流切应力的作用，泥沙颗粒不断起动，在湍动扩散的作用下，泥沙在不断向前运动的同时，也向上扩散，〈s〉/sb沿程逐渐增大.

图8中给出了算例2长槽道中瞬时和时均的断面平均泥沙浓度的沿程分布，从图中可以看出，虽然经过断面平均，瞬时的泥沙浓度仍然有较大的脉动.但在总体趋势上，瞬时的断面平均泥沙浓度呈现与时均的断面平均的泥沙浓度一致的趋势，泥沙浓度随着x的增大而逐渐变大.从图中可以看出，从x=0H开始，断面平均泥沙浓度〈s〉/sb沿着x方向逐渐增大，但增大速度逐渐减慢，在x=42H处，达到最大值0.14，但由于断面平均泥沙浓度〈s〉/sb关于x的导数大于零，整个长槽道中，床面一直处于净冲刷状态.

图8 算例2中断面平均的泥沙浓度沿程分布Fig.8 Distribution of the cross-section averaged sediment concentration in case 2

泥沙浓度湍动强度沿水深的分布见图9(a).图中横坐标z采用摩阻流速u∗和黏性系数υ进行无量纲化，有z+=zu∗/υ，采用平衡时近底部的泥沙浓度sb将泥沙湍动强度s′进行无量纲化.从图中可以看出泥沙的湍动强度s′/sb在近底面处达到最大值0.49，之后随着水深的增大而急剧降低，在水面附近仅为0.04，在自由水面处几乎为零.图9(b)给出了采用LES得到的垂向上泥沙的湍动通量沿水深的分布.横坐标为z+=zu∗/υ，纵坐标w′s′采用u∗sb进行无量纲化，在近底部w′s′/u∗sb取得最大值0.115，之后随着水深的增大不断减小，在水面处基本为零.

图9 算例1中泥沙浓度湍动强度和湍动通量的分布Fig.9 The turbulence intensity of sediment concentration and vertical sediment turbulence flu in case 1

5.3 湍动黏性系数与扩散系数分布

对于湍动切应力和通量有

其中vt，εt分别为湍动黏性系数和湍动扩散系数，是湍动切应力是垂向湍动泥沙通量.

图10中给出了循环槽道中(算例1)和长槽道中(算例2)采用LES得到的湍动黏性系数和扩散系数沿着水深的分布.图中的横坐标z采用水深H进行无量纲化，纵坐标湍动黏性系数υt和湍动扩散系数 εt，分别采用分子扩散系数 υ和 ε进行无量纲化，υ=1/20000，取分子Schmidt数为1，因此ε=1/20000.湍动扩散系数εt理论值采用式(9)计算

其中κ为卡帕常数，u∗为摩阻流速，z为离开床面的距离，H为水深.

从图10(a)中可以看出，采用大涡模拟计算得到的湍动黏性系数υt/υ关于z/H=0.5不对称分布，在z/H=0.35处，υt/υ取得最大值60，而大涡模拟计算得到的扩散系数εt/ε关于z/H=0.5对称分布，在水深中间处，εt/ε取得最大值70.在z/H介于[0.15, 0.85]的区间中时，εt/ε的理论值明显大于LES计算值.在水深中间处由理论公式得到的εt/ε的最大值是LES计算得到的εt/ε最大值的1.4倍.这是因为理论公式推导时假定：(1)流速服从对数分布；(2)泥沙湍动扩散系数等于动量扩散系数.而实际上流速只在对数区服从对数分布，由于水流中泥沙的跟随性不如普通的污染物，计算得到的湍动扩散系数与理论值有差别.

图10 湍动黏性系数和扩散系数分布图Fig.10 The turbulence viscosity coefficient and the turbulence dif f usion coefficient in two cases

图10(b)给出了长槽道中(算例2)x=10H,20H, 30H,40H处展向平均的湍动黏性系数和扩散系数的分布.从图中可以看出，由于统计样本数较少，因此动黏性系数υt和扩散系数εt有轻微的波动，但总体上看，黏性系数υt和扩散系数εt沿着x向的各个断面基本相同，并且沿着垂向呈现先变大后变小的趋势.

图11 算例1中湍动Schmidt数分布图Fig.11 The profil of the turbulence Schmidt number in case 1

图11给出了湍动 Schmidt数 σt沿着水深的分布.从图中可以看出σt沿着水深不断变化，湍动Schmidt数σt在接近水面和底部时较大，在水深中间处较小，当z/H介于[0.25,0.85]中时，σt小于1.在z/H=0.55时，σt=0.83为最小值.但在雷诺时均模拟中σt常常假定为常数[30-31]，当水流结果计算准确时，σt存在偏差，从而在泥沙浓度计算中引入误差，尤其对近底部和水面附近的泥沙浓度有较大的影响.

6 结论

在本文中采用基于并行计算的大涡模拟程序，借鉴了雷诺时均模拟中悬沙计算中的参考高度概念，建立了LES中的细颗粒泥沙净冲刷条件下的近底边界条件模拟了循环槽道和长槽道中水流运动和泥沙输移的过程.在输沙平衡条件下，泥沙的浓度分布服从Rouse公式的分布.长槽道流动中，流动沿着x向均匀分布，整个槽道处于净冲刷状态，在水流切应力和湍动作用下，泥沙颗粒不断起动并向上扩散，泥沙浓度沿程不断增大.泥沙浓度的湍动强度和垂向湍动通量在近底处达到最大，沿着垂向迅速衰减.通过LES统计得到的湍动黏性系数和扩散系数沿着水深方向呈现先增加后减小的趋势，整体上湍动扩散系数的计算值小于实际值.湍动Schmidt数不是常数，当z/H介于[0.25，0.85]中时，σt小于1.在z/H=0.55时，σt=0.83为最小值.

1 Walling DE,Owens PN,Carton J,et al.Storage of sedimentassociated nutrients and contaminants in river channel and flood plain systems.Applied Geochemistry,2003,18(1)：195-220

2 Collins AL,Walling DE,Leeks GJ.Storage of fine-grainesediment and associated contaminants within the channels of lowland permeable catchments in the UK//Sediment Budgets 1,Walling D E,Horowitz A J,eds.IAHS Publication No.291,IAHS Press：Wallingford,UK,2005：259-268

3 Fang HW,Zhao HM,Shang QQ,et al.Ef f ect of biofil on the rheological properties of cohesive sediment.Hydrobiologia,2012, 694(1)：171-181

4 Huang L,Fang H,Reible D.Mathematical model for interactions and transport of phosphorus and sediment in the Three Gorges Reservoir.Water Research,2015,85：393-403

5 Fang HW,Huang L,Wang JY,et al.Environmental assessment of heavy metal transport and transformation in the Hangzhou Bay, China.Journal of Hazardous Materials,2015,302(17)：3836-3849

6 Marchioli C,Soldati A.Mechanisms for particle transfer and segregation in a turbulent boundary layer.Journal of Fluid Mechanics, 2002,468：283-315

7 Chang YS,Scotti A.Entrainment and suspension of sediments into a turbulent fl w over ripples.Journal of turbulence,2003,4(19)：1-22

8 van Rijn LC.2007.Unifieview of sediment transport by currents and waves II：suspended transport.Journal of Hydraulic Engineering-ASCE,2007,133(6)：668-689

9 Zeng J,Constantinescu G,Weber L.A 3D non-hydrostatic model to predict fl w and sediment transport in loose-bed channel bends.Journal of Hydraulic Research,2008：46(3)：356-372

10 Fang H,He G,Liu J,et al.3D numerical investigation of distorted scale in hydraulic physical model experiments.Journal of Coastal Research,2015,10052(1)：41-54

11 Keylock CJ,Hardy RJ,Parsons DR,et al.The theoretical foundations and potential for large-eddy simulation(LES)in fluvia geomorphic and sedimentological research.Earth-Science Reviews, 2005,71(3-4)：271-304

12 Rodi W.Comparison of LES and RANS calculations of the fl w around bluf fbodies.Journal of Wind Engineering and Industrial Aerodynamics,1997,69-71：55-75

13 Zedler EA,Street RL.Large-eddy simulation of sediment transport：currents over ripples.Journal of Hydraulic Engineering-ASCE, 2001,127(6)：444-452

14 Zedler EA,Street RL.Sediment transport over ripples in oscillatory fl w.Journal of Hydraulic Engineering-ASCE,2006,132(2)：180-193

15 Chou YJ,Fringer OB.Modeling dilute sediment suspension using large-eddy simulation with a dynamic mixed model.Physics of Fluids,2008,20(11)：115103

16 Widera P,Ghorbaniasl G,Lacor C.Study of the sediment transport over fla and wavy bottom using large-eddy simulation.Journal of Turbulence,2009,10(33)：1-20

17 Widera P,Toorman E,Lacor C.Large eddy simulation of sediment transport in open-channel fl w.Journal of Hydraulic Research, 2009,47(3)：291-298

18 Bai J,Fang HW,Stoesser T.Transport and deposition of fin sediment in open channels with dif f erent aspect ratios.Earth Surface Processes and Landforms,2013,38(6)：591-600

19 Germano M,Piomelli U,Moin P,et al.A dynamic subgrid-scale eddy viscosity model.Physics of Fluids,1991,3(7)：1760-1765

20 Breuer M,Rodi W.Large-eddy simulation of turbulent fl w through a straight square duct and a 180°bend//Voke PR,Kleiser L,Chollet JP,eds.Direct and Large-Eddy Simulation I.Netherlands：Kluwer, 1994.273-285

21 Hinterberger C,Frohlich J,Rodi W.Three-dimensional and depthaveraged large-eddy simulations of some shallow water fl ws.Journal of Hydraulic Engineering-ASCE,2007,133(8)：857-872

22 Stoesser T,Nikora V.Flow structure over square bars at intermediate submergence：large eddy simulation(LES)study of bar spacing ef f ect.Acta Geophysica,2008,56(3)：876-893

23 Stoesser T,Ruether N,Olsen NRB.Calculation of primary and secondary fl w and boundary shear stresses in a meandering channel.Advances in Water Resources,2010,33(2)：158-170

24 Zhu J.A low-dif f usive and oscillation-free convection scheme.Communications in Applied Numerical Methods,1991,7(3)：225-232

25 Stone HL.Iterative solution of implicit approximations of multidimensional partial dif f erential equations.SIAM Journal on Numerical Analysis,1968,5(3)：530-558

26 CelikI,RodiW.Modelingsuspendedsedimenttransportinnonequilibrium situations.Journal of Hydraulic Engineering-ASCE,1988, 114(10)：1157-1191

27 van Rijn LC.Sediment transport,part II：suspended load transport.Journal of Hydraulic Engineering-ASCE,1984,110(11)：1613-1641

28 Werner H,Wengle H.Large-eddy simulation of turbulent fl w over and around a cube in a plate channel//Proceeding of 8th Symposium on Turbulent Shear Flows,1991.155-168

29 Mito Y,Hanratty TJ,Zandonade P,et al.Flow visualization of super bursts and of the Log-Layer in a DNS at Reτ=950.Flow Turbulence and Combustion,2007,79(2)：175-189

30 Wu W,Rodi W,Wenka T.3D numerical modeling of fl w and sediment transport in open channels.Journal of Hydraulic Engineering-ASCE,2000,126(1)：4-15

31 Fang HW,Wang GQ.Three-dimensional mathematical model of suspended-sediment transport.Journal of Hydraulic Engineering-ASCE,2000,126(8)：578-592

NUMERICAL SIMULATION OF EROSION AND TRANSPORT OF FINE SEDIMENTS BY LARGE EDDY SIMULATION

Bai Jing∗Fang Hongwei†,1)He Guojian†Xie Chongbao∗Gao Hong∗∗

∗(China Irrigation and Drainage Development Center，Beijing100054，China)

†(State Key Laboratory of Hydro Science and Engineering,Department of Hydraulic Engineering，Tsinghua University，Beijing100084，China)

∗∗(China Green Water International Consulting Co.，Ltd，Beijing 100054，China)

In general Reynolds-averaged simulation(RANS)is used in the traditional numerical simulation of water fl w and sediment transport.Large eddy simulation(LES)can reflec fl w structures more accurately and give more details of water fl w compared with RANS.The development of computing technology makes it possible to study the rules of water fl w and sediment transport by an LES model.In this paper,we tried to introduce boundary conditions for suspended sediment transport for the LES model under the net-erosion condition.Water fl w and sediment transport in a cyclic case and a one-way fl w case were calculated via the LES model with a dynamic sub-grid stresses module and a suspended sediment calculation module in the paper.Direct numerical simulation(DNS)results were used to calibrate the LES model and the results from LES showed good agreements with the DNS results.The distribution characteristics of sediment concentration,turbulence intensity and turbulent flu es of sediment were explored in the paper.Under the net-erosion condition,the equilibrium sediment concentration profil was coincident with the line of the Rouse equation.It showed that the turbulence intensity and turbulent flu es of sediment had peak values near the bottom and then decreased rapidly along the vertical direction.The turbulent viscosity and dif f usion coefficients were calculated and their peak values were or near the mid-depth of water.The turbulent Schmidt number was not constant along the vertical direction,and it was larger near the free surface and the bottom while it was smaller near the mid-depth of water fl w.

fin sediment,net erosion,LES,sediment transport

TV142

A doi：10.6052/0459-1879-16-235

2016-08-25收稿，2016-11-15录用,2016-11-18网络版发表.

1)方红卫，教授，主要研究方向：环境泥沙及河流动力学.E-mail：fanghw@mail.tsinghua.com

白静，方红卫，何国建，谢崇宝，高虹.细颗粒泥沙净冲刷和输移的大涡模拟研究.力学学报，2017,49(1)：65-74

Bai Jing,Fang Hongwei,He Guojian,Xie Chongbao,Gao Hong.Numerical simulation of erosion and transport of fin sediments by large eddy simulation.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：65-74