孤立波作用下边界层的模拟与分析

张 健，刘昭伟，包洪福

（1.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；

2.清华大学 水利水电工程系，北京 100084；3.中国长江三峡集团公司，北京 100038）

孤立波作用下边界层的模拟与分析

张 健1，2，刘昭伟1，2，包洪福2，3

（1.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；

2.清华大学 水利水电工程系，北京 100084；3.中国长江三峡集团公司，北京 100038）

孤立波作用下的边界层内的剪切力以及涡量变化对海啸传播和海底地形塑造十分重要。本文基于多区域谱方法，利用直接模拟（DNS）数值模型，对在具有矩形断面的U形水洞内的孤立波下的边界层流动进行了模拟。将数值模拟结果与解析解以及试验结果进行了对比，发现数值结果与后两者吻合得较好。模拟结果显示，在低雷诺数下，扰动不会改变流态，而随着雷诺数的增大，流态会变得十分复杂。中等雷诺数情况下，边界层内会产生正向的涡，并进行稳定的传播。在较高雷诺数情况下，流动进入层流向紊流发展的过渡期，此时边界层内会产生正涡以及负涡，并会在水深方向进行不规则运动。

孤立波；边界层；剪切力；涡管

1 研究背景

海啸波进入大陆架浅水区时，会受到水深急剧变化的影响，从而能量瞬时集中，波高骤然增大。在这一过程中非线性作用突显，传统线性Stokes波浪理论已不再适用，需要采用孤立波模型对其进行研究。在高雷诺数假设下，远离水底的水体内黏性作用很小，可忽略不计，因而可以利用有势流动的解析方法进行分析。但海啸波的长时间、长距离传输过程则会对其进行精确预测产生非常明显的负面作用，如传播过程中的阻力影响难以估算、边界层分离影响程度等等。因而需要利用孤立波模型对这些问题进行进一步的探究。

通过假设孤立波作用下边界层内流动为层流状态，可以对这一问题进行理论分析。例如，Liu and Orfila利用势流理论推导了瞬时长波传播时，边界层内的流速分布以及剪切力分布的解析解。［1］Liu等通过迭代的方法对边界层内非线性的流动方程进行求解，发现非线性方程结果与线性方程结果相差很小，得到了在一定情况下可以只考虑线性项在边界层内的影响的结论。［2］

除理论分析外，数值模拟也是对边界层进行探究的有效手段。Vittoriand Blondeaux利用RANS方法对在自由表面的孤立波作用下的底部边界层进行模拟。结果显示在减速期、波高大于某一临界值情况下会出现紊动现象，该临界值与边界层厚度与当地水深之比有关，而且边界层内剪切系数与多种参数有关［3］。Suntoyo and Tanaka同样利用RANS方法对孤立波边界层问题进行探究，采用BSL k-ω紊流模型计算线性化后的边界层方程。模拟结果显示随着雷诺数的增加，紊动能以及底面剪切力会随之增加，但反向剪切力的峰值却会减小。在完全发展的紊流状态下，边界层内剪切力与边界层外紊动扩散所产生的速度同时产生，而且反向剪切力峰值同样减小，这一结论对近岸泥沙输移影响巨大［4］。

利用试验可以对理论分析和数值计算进行验证。Liu等在波浪水槽中进行了一系列孤立波下边界层速度场测量试验，试验中孤立波由活塞式造波机产生，并利用PIV对速度场进行测量。试验条件下边界层为层流状态。试验结果与Liu and Orfila和Liu etal对层流流场的分析结论相一致。［2］

仅利用活塞制造的波动很难在底部形成紊动，为对边界层内的相干结构进行探究，Sumer et al在具有矩形断面的U形水洞（water tunnel）中进行了一系列探究性的试验。试验利用气压系统模拟制造孤立波运动形态，利用水洞是U形的这一特殊构造，模拟孤立波中流体微团的速度变化，从而使底部出现不稳定状态，边界层内不再保持层流［5］。试验得到了雷诺数介于9.4×104和2.0×106之间的工况下的边界层内速度和底面剪切力的数据。此处雷诺数的定义为：

式中：U0m为边界层外自由流动的最大速度，m/s；ν为运动黏滞系数，m2/s；2a为孤立波作用下自由流动区域中一个周期相位内流体质点的最大位移的1/2，m。

Sumer et al对试验现象进行分析，发现当 Re≤2×105时，流场为层流状态；当 2×105＜Re≤5×105时，流动基本保持层流状态，但当压力梯度变向时会在边界层内出现一系列涡管，涡管的出现，使得底面剪切力分布也变得不再光滑；当雷诺数达到5×105时，流场由层流开始向紊流转变，进入过渡期。当雷诺数达到试验条件下的最大值2.0×106情况时，约有半个周期进入了紊流状态，而开始的半个周期则仍保持层流状态。根据雷诺数增长引起的边界层内剪切力变化规律来看，雷诺数继续增加会导致紊动相位持续的时间持续增加，甚至遍布几乎整个周期，即流动进入完整的紊动流动状态，但现在试验还无法确定这一位置的具体雷诺数情况。［5］

有别于Vittori和 Blondeaux以及Suntoyo和 Tanaka的大涡模拟方法，本文采用DNS法，利用Di⁃amessis etal提出的多区域谱方法模型［6］，对Liu和Orfila、Liu et al的理论分析以及Sumer et al的试验分析进行验证，并对不同时期的边界层内流场分布的特性进行探究。

2 模型建立

2.1 控制方程二维U型水槽中的流动可以用Navier-Stokes方程进行描述：

其中，u=（u，w）为速度矢量；p为压力项，Pa；ρ为水的密度，kg/m3，假设水体不可压缩，则ρ可以取常数1.0×103kg/m3；ν为运动黏滞系数，取值为1.0×10-6m2/s。

在边界层内由于黏性作用而无法保持有势运动，根据这一特性，我们可以将边界层内的压强场分解为两部分进行求解，即有势流动引起的平均项和有旋运动引起的扰动项，这样压强场可以表示为：

其中下标为w的表示的是平均项；带有上标波浪线的为扰动项。压强平均项利用Liu和Orfila的解析方法进行推导

2.2 失稳触发通过对非恒定流的振荡边界层进行测量以及数值模拟显示紊动产生的动力学条件十分复杂。Costamagna，Vittor和Blondeaux认为在振荡流动的水体中产生紊动现象与稳定水体中的紊动现象有相同的机理［7］。在加速期末段的近壁区域中会出现低速带，而在减速期这些低速带就会扭转、振荡，最终破碎，产生小型涡。由于黏性作用，这些涡又会在加速期的初段耗尽。因而一般的紊流模型只能显示不稳定边界层的整体特性，却不能显示这一过程的变化特征。

在Sumer etal的U型水洞试验中，我们可以看到边界层内的涡流管会沿一定距离分布，并同时产生和消失，换言之涡流管的生成和发展具有一定周期性。在数值模拟中，如果没有引入任何扰动的话，流场几乎不会进入失稳状态。为了使数值模拟中的流场迅速进入失稳状态，并可以对不稳定的变化过程有更清晰的探究，可以在最易失稳的相位下引入一个微小的速度扰动分量，此相位下流动速度达到最大值，压力梯度开始变向。速度扰动分量的纵向和垂向表达形式为：

其中，m表示扰动速度的峰值，m值要大到足够引起流场失稳，也要小到不会继续对失稳后的流场产生影响，因而在本文的模拟中m=0.001；b表示初阶不稳定的相对强度，Smyth和 Moum认为b=0.4177［8］；h0表示边界层的厚度；z0表示在0°相位位置时边界层的垂向坐标；u0=U0m；k0为沿纵向方向的波数。

本文所采用的数值模型在纵向方向上采用切比雪夫多项式进行拟合，因而扰动速度最终表达形式也是以三角级数的形式展现出来，而且可以式（6）、式（7）满足连续方程以及下文的边界条件。

为了计算h0、z0以及k0三个参数，需要引入剪切层涡厚度这一参数，其定义为：

δω虽然会随相位改变而有所不同，但差异并不大，可以用ωt≈32°的位置计算涡厚度。于是上述参数的计算公式为：

式中λ表示波长。

2.3 边界条件与初始条件根据Sumer et al的试验，我们将计算域定义在一个Lx×H的矩形区域内，其中H为U型水槽从槽底到中心线的高度；而Lx为沿波流动方向的可变长度，对于不同试验有所不同，主要由涡流管的个数和尺寸决定，以便满足速度和压力的周期性特性：

计算域的顶部为无滑移、无变形的边界：

底部为固定边界：

2.4 数值方法控制方程利用Diamessis et al提出的多区域谱方法模型进行求解［6］。模型中时间导数项利用向后分步差分进行离散，共分3步交替进行，以保证时间计算的准确性。二维模型中两个方向具有不同的性质，因而也采用不同的离散方法对空间导数项进行离散。在周期性方向（波运动方向）采用傅立叶谱法进行离散，基函数选用切比雪夫多项式，模态阶数根据不同问题进行选取（本文计算案例中采用前8阶）。在非周期性方向（水深方向），将计算域划分为多个子区域，每个区域的垂向长度不同，为更好的对边界层进行模拟，由下至上区域高度逐渐提高。计算时，每个子区域内采用勒让得谱法进行离散，模态阶数同样由问题精度要求决定（本文计算案例中采用了前4阶）。较高雷诺数以及较薄边界层的情况下，子区域间数值传递往往会出现Gibbs振荡，为避免误差的继续发展，模型依次采用了罚函数、谱过滤以及界面平均三种方法来保证计算的稳定以及误差最小化。

本文的模拟都从静止状态出发：

3 算例分析

本文共对4种工况进行模拟，分别对应三种不同的流场状态：无涡层流、有涡层流和过渡期。利用前文所建立的数值模型，对不同流态的底面剪切力、涡量分布等进行分析，并与理论分析和试验结果进行比较。表1为4种工况的参数设置情况，各个参数的物理含义如前文所述。

表1 算例参数

对于具体的算例来说，纵向计算长度为Lx，垂向长度为0.145m。计算时为了能更好地反映边界层以及水体底部的信息，在垂向上将计算域划分为11个子区域，由下至上高度分别为0.01、0.02、0.04、0.08、0.16、、0.32、0.64、1.28、2.56、5.12以及4.27 cm。经过网格无关性分析之后发现纵向96个节点，垂向每个子区域内25个节点就可以对问题进行较为精确地模拟，这种划分方式可以使得效率和精度得到较好平衡。

3.1 无涡层流式（6）、式（7）表示的扰动速度对涡管的迅速产生有重要的影响。但在低雷诺数情况下，引入扰动是否也会产生涡管，则需要进一步的分析。本节将对表1的工况1进行模拟来说明无涡层流的流场分布情况。数值模拟时网格参数如上文所述。

低雷诺数情况下底面剪切力分布如图1所示。图中实线部分为数值计算结果，虚线部分为Liu和 Orfila、Liu等提出的解析解，两者对比可以发现具有较好的吻合度。而且在任何相位位置都未发生震荡现象，流动保持稳定，直至速度降到0。

图1 Re=2.56×104时底面剪切力数值结果与解析结果相位对比

图2为数值计算得到的3种不同相位下涡量分布等势图。从图中可以看出流场未出现涡管，维持层流状态，与剪切力分布相互验证。而且，随着相位的发展，涡量也逐渐减小，直至变为0，满足层流特性。

3.2 有涡层流Sumer等根据试验结果将孤立波作用下边界层内的流动分为4个阶段，并给出了每个阶段具体的雷诺数分布情况［5］，但实际上两个阶段之间很难有较为明显的界限，尤其是有涡管层流与过渡区之间的界限更加不明显。利用工况2、工况3进行对比来分析中等雷诺数情况下边界层内流动情况。与3.1节不同的是，进入不稳定阶段后很难进行理论求解，因而将数值计算结果与试验进行对比。

图3分别为Re=4.20×105，Re=6.06×105时底面剪切力分布数值解与试验结果对比情况。从图中可以看出，中等雷诺数下在加速期过后流动进入不稳定状态，剪切力出现轻微震荡。由于这种震荡具有随机性，因而试验结果与数值解果略有差距，但二者的趋势却是一致的。另一方面，有涡层流和过渡期的界限并不十分明显，但雷诺数越大，震荡幅度越大，延续的相位也越大。

图2 Re=2.56×104时3种相位时边界层内涡量分布等势图

图3 中等雷诺数时底面剪切力数值解与试验点相位对比

图4 Re=4.20×105时3种相位时边界层内涡管分布等势图

图5 Re=6.06×105时3种相位时边界层内涡管分布等势图

图4、图5分别为Re=4.20×105，Re=6.06×105时边界层内涡管发展等势图。从图中可以看出，中等雷诺数情况下，涡管产生后会与波进行反向传播，但涡管中心高度不会改变，并稳定的在边界层内发展。而有涡层流与过渡前期的涡管分布区别并不大。

3.3 过渡期随着雷诺数进一步的增加，边界层内流动会发生巨大的变化。接下来将对层流向紊流转变的过渡期进行模拟分析，所选工况为工况4。

图6为大雷诺数下底面剪切力分布情况。与中等雷诺数情况对比可以看到，此时剪切力的震荡情况更加剧烈，且振幅已经超过了稳定状态下由于压力梯度引起的最大剪切力。试验结果与数值模拟结果有相同的趋势，但相位略有差距。随着雷诺数的增加，试验测量的难度也急剧增大，因而这种情况下的测量结果也很难与实际情况完全相同，因而可以认为数值结果在一定程度上是合理的。

图6 Re=1.22×106时底面剪切力数值解与试验点相位对比

图7为工况的涡管分布情况。从图中可以看出，大雷诺数下涡管的运动趋势与低雷诺数相似，但复杂程度却加大了，涡管的分布有逐渐上移的趋势，这使得边界层的界限也跟着上移。

图7 Re=1.22×106时3种相位时边界层内涡管分布等势图

4 结论与讨论

本文利用Diamessis etal.提出的多区域谱方法模型，对孤立波作用下边界层内流动进行了分析，并将计算结果与Liu和Orfila、Liu et al的解析解以及Sumer etal的试验分析进行了对比。通过数值模拟得到：（1）利用二维模型可以对U型水槽试验进行模拟，模拟结果也较为理想；（2）低雷诺数（Re≤2×105）情况下流动不会受到扰动的影响，可以完全保持层流状态，边界层内并不产生涡管；（3）中等雷诺数（2×105＜Re≤5×105）情况下，边界层内会产生涡管，但流动仍可以维持在层流状态，底面剪切力会小幅振荡，但会随相位逐渐减小，涡产生后会与波进行反向运动，但涡管中心一般不会改变，流动较为稳定；（4）高雷诺数（Re＞5×105）情况下，流动进入过渡状态，底面剪切力发生剧烈振荡，振幅较大，涡管的产生与发展出现复杂情况。

自然状态下，流动多处于高雷诺数情况，而低、中雷诺数则只在波产生的初期出现，因而关于高雷诺数情况下涡管产生、发展以及相互作用仍需要做进一步的探究。

参 考 文 献：

［1］ Liu P L-F，Orfila A.Viscous effects on transient long-wave propagation［J］.J.Fluid Mech.，2004，520：83-92.

［2］ Liu P L-F，Park Y S，Cowen E A.Boundary layer flow and bed shear stress under a solitary wave［J］.J.Fluid Mech.，2007，574：449-463.

［3］ Blonderaux P，VittoriG.RANSmodelling of the turbulent boundary layer under a solitary wave［J］.Coastal Engi⁃neering，2012，60：1-10.

［4］ Suntoyo，Tanaka H.Numericalmodeling of boundary layer flow for a solitarywave［J］.J.Hydro-environmentRe⁃search，2009，3：129-137.

［5］ Sumer B M，Jensen PM，Sorensen L B，et al.Coherent structures in wave boundary layers.Part2.solitarymo⁃tion［J］.J.Fluid Mech.，2010，646：207-231.

［6］ Diamessis P J，Domaradzki JA，Hesthaven JS.A spectralmultidomain penaltymethod model for the simulation of high Reynolds number localized incompressible stratified turbulence［J］.J.Compu.Phys.，2004，202（1）：298-322.

［7］ Costamagna P，VittoriG，Blondeaux P.Coherent structures in oscillatory boundary layers［J］，J.Fluid Mech.，2003，474：1-33.

［8］ Smyth W D，Moum JN.Length scales of turbulence in stably stratified mixing layers［J］.J.Phys.Fluids.，2000，12（6）：1327-1342.

Simulation and analysis of boundary layer under the solitary wave

ZHANG Jian1，2，LIU Zhaowei1，2，BAO Hongfu2，3
（1.State Key Laboratory of Hydrology-Water Resourcesand Hydraulic Engineering，HohaiUniversity，Nanjing 210098，China；2.Departmentof Hydraulic Engineering，Tsinghua University，Beijing 100084，China；3.China Three GorgesCorporation，Beijing 100038，China）

The shear stress and vortex tube in the boundary layer are very important to tsunami propaga⁃tion and underwater terrain under the solitary wave.This paper takes advantage of direct numerical simula⁃tion method（DNS） to simulate the boundary layer under the solitary wave in a U shape water tunnel，based on the spectral multi-domain model.By comparing numerical results with analytical solutions and ex⁃perimental results，it is found that they agree with each other well.The numerical results show that with a low Reynolds number，the disturbance does not change the flow pattern，and with the increase of Reynolds number，flow pattern will become very complicated.Within a medium Reynolds number，positive vortex tube will appear in the boundary layer，and transport stably.Under the condition of a high Reynolds num⁃ber，layer flow changes into a transition period of turbulence flow and positive vorticity and negative vortici⁃ty appear in the boundary layer，and have oscillating development.

solitary wave；boundary layer；shear stress；vortex tube

TV139.2

：Adoi：10.13244/j.cnki.jiwhr.2015.04.003

1672-3031（2015）04-0254-07

（责任编辑：李福田）

2015-04-30

国家自然科学基金资助项目（51279079）；“十二五”国家科技支撑计划（2013BAB05B05，2010BAE00739）

张健（1993-），男，吉林人，硕士生，主要从事水力学研究。zhang0249@126.com

刘昭伟（1973-），男，河北人，副教授。liuzhw@tsinghua.edu.cn