垂向基于谱方法的三维弯道水流模型
杨飞, 傅旭东     
清华大学 水利水电工程系, 北京 100084
摘要:传统平面二维水动力模型模拟弯道水流时,将流速垂向分布固定化,不能随边界变化和二次流发展进行调整,而垂向完全自动调整的三维水动力模型计算效率低。该文采用谱方法将平面流速的垂向分布用正交多项式表示,通过加权余量法建立简易三维模型,对流项采用基于Gauss点的半Lagrangian法(semi-Lagrangian scheme,SLS)计算。由急弯弯道水流试验验证可知,多项式阶数大于1时模拟结果较好。弯道主流线位置预测整体误差不超过7%,精度和传统的三维水动力学模型相当。对涡黏性系数简单求解,但考虑了垂向和水平方向的差异性,能够得到合理的弯道水流结构。垂向基于谱方法建立的三维模型,无垂向网格离散,计算量明显减少,获得了与平面二维模型相当的计算效率。
关键词谱方法    三维水动力学模型    二次流    弯道水流    各向异性    
3-D hydrodynamic model using the spectral method in the vertical direction for bend flow simulations
YANG Fei, FU Xudong     
Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China
Abstract: The traditional depth-averaged 2-D hydrodynamic models for bend flow simulations assume velocity profiles for the secondary flow terms in the momentum equations which are unable to adjust to dynamic conditions. More flexible three-dimensional models are not very efficient. A simplified 3-D model using a spectral method in the vertical direction was developed with the flow velocity components modeled by orthogonal polynomials in the vertical direction using polynomial coefficient equations obtained using the weighted residuals method with the advection terms defined at the vertical Gauss points by the semi-Lagrangian scheme. Simulated flow structures in a sharp bend open channel match well with measured data for polynomials having degrees larger than 1 with reasonable flow structures. The mean error of the predicted main flow location is less than 7%, equivalent to other 3-D hydrodynamic models. The eddy viscosity is solved in a simple way with consideration of the turbulence anisotropy between the vertical and horizontal directions. Since this method does not have a vertical grid, the calculational efficiency is proved to be to 2-D models.
Key words: spectral method     3-D hydrodynamic model     secondary flow     bend flow     anisotropy    

常见的冲积河流河型中的水流流路多发生不同程度的弯曲,弯曲流动的离心力和横向水面比降不匹配产生弯道二次流(螺旋流),它控制了河流横向的动量、泥沙输运和长期的横向演变趋势。对于弯曲河流模拟,三维模型不需要假定垂向水流结构直接解析二次流,但计算效率低;垂向平均的二维模型计算效率高,能满足一般河流水沙和河床冲淤计算的精度要求,是常用的模型工具。平面二维模型模拟弯曲河流流速、床沙、冲淤等特征的横向分布时需要考虑二次流修正。

一类二次流修正方法[1-2]都是在微弯R/B < 0.5(R为曲率半径,B为河宽)情况下,不考虑二次流对主流的反馈作用,称为线性方法。急弯(R/B>0.5)情况下,二次流对主流的影响加强,二次流强度不再线性增加[3-5],需考虑反馈作用,称为非线性方法[5-11]。另外,又有考虑水流结构调整滞后于弯曲边界的变化的二次流松弛模型[12-14]

无论线性方法还是非线性方法,都是将参数化垂向水流结构作为额外的弥散应力项(dispersion stresses)[4]添加到垂向积分的动量方程中,虽然简化了二次流计算,但降低了模型的自由度,使得模型对边界条件以及流速的垂向分布的调整变化不能完全响应。因此如何使模型能够像三维模型那样计算流速的垂向分布,实时反映流场变化,同时避免三维模型的复杂计算是本文的研究内容。在此引入谱方法的思想,将流速的垂向分布采用多项式表示,实时计算多项式的系数调整。高效的谱方法对光滑函数指数型逼近有谱精度,能够得到高精度的数值解。目前谱方法在计算流体动力学中得到较大的发展[15],但还没有应用到三维浅水方程的计算中。

为了弥补平面二维模型二次流方法适应性差的缺陷,同时避免传统三维模型垂向离散带来的大量计算,本文给出了一个基于平面二维网格与垂向低阶谱方法结合的模型。用低阶多项式拟合流速的垂向分布,避免了垂向离散,求解加权余量法得到的变分方程,实现流速的垂向分布动态地响应边界条件变化。水槽试验验证表明,采用简单的紊流模型通过考虑涡黏性系数的各向异性,弯道主流线和水流流速分布模拟结果与实测一致。结合计算效率与精度,该方法有很强的实际意义。

1 模型建立 1.1 控制方程

控制方程采用曲线坐标系下的三维浅水方程。曲线坐标系采用平面二维曲坐标和垂向z坐标,并假定垂向与平面正交。对于一般曲线坐标系,可以参考文[16]。

连续性方程

$ \frac{\partial }{{\partial \xi }}\frac{u}{J} + \frac{\partial }{{\partial \eta }}\frac{v}{J} + \frac{\partial }{{\partial z}}\frac{w}{J} = 0. $ (1)

水平动量方程

$ \begin{array}{l} \frac{{\partial u}}{{\partial t}} + J\left( {\frac{\partial }{{\partial \xi }}\frac{{uu}}{J} + \frac{\partial }{{\partial \eta }}\frac{{uv}}{J} + \frac{\partial }{{\partial z}}\frac{{uw}}{J}} \right) = \\ - {\alpha _1}uu - {\alpha _2}uv - {\alpha _3}vv - \\ g\left( {{\beta _1}\frac{{\partial H}}{{\partial \xi }} + {\beta _2}\frac{{\partial H}}{{\partial \eta }}} \right) + \frac{\partial }{{\partial z}}\left( {{v_{\rm{z}}}\frac{{\partial u}}{{\partial z}}} \right) + \\ \frac{\partial }{{\partial \xi }}\left( {{v_{\rm{h}}}\xi _r^2\frac{{\partial u}}{{\partial \xi }}} \right) + \frac{\partial }{{\partial \eta }}\left( {{v_{\rm{h}}}\eta _r^2\frac{{\partial u}}{{\partial \eta }}} \right), \end{array} $ (2a)
$ \begin{array}{l} \frac{{\partial v}}{{\partial t}} + J\left( {\frac{\partial }{{\partial \xi }}\frac{{uv}}{J} + \frac{\partial }{{\partial \eta }}\frac{{vv}}{J} + \frac{\partial }{{\partial z}}\frac{{vw}}{J}} \right) = \\ - {\alpha _4}uu - {\alpha _5}uv - {\alpha _6}vv - \\ g\left( {{\beta _3}\frac{{\partial H}}{{\partial \xi }} + {\beta _4}\frac{{\partial H}}{{\partial \eta }}} \right) + \frac{\partial }{{\partial z}}\left( {{v_{\rm{z}}}\frac{{\partial v}}{{\partial z}}} \right) + \\ \frac{\partial }{{\partial \xi }}\left( {{v_{\rm{h}}}\xi _r^2\frac{{\partial v}}{{\partial \xi }}} \right) + \frac{\partial }{{\partial \eta }}\left( {{v_{\rm{h}}}\eta _r^2\frac{{\partial v}}{{\partial \eta }}} \right). \end{array} $ (2b)

其中:t为时间;H为水位;h为水深;g为重力加速度;xyz为直角坐标;ξη为平面曲坐标;νzνh为垂向与水平涡黏性系数;J=ξxηy-ξyηx为Jacobian行列式;uv为水平流速在曲坐标系下的协变分量;w为垂向流速;ξrηr为曲坐标系两个方向Lame系数的倒数;ζ=2(z-zb)/h-1为相对高度;zb为床面高程。动量方程中的水平扩散项表达式已做近似简化处理,压力梯度项系数

$ \begin{array}{l} {\beta _1} = \xi _x^2 + \xi _y^2,{\beta _2} = {\xi _x}{\eta _x} + {\xi _y}{\eta _y},\\ {\beta _3} = {\xi _x}{\eta _x} + {\xi _y}{\eta _y},{\beta _4} = \eta _x^2 + \eta _y^2. \end{array} $ (3a)

对流项坐标转换产生项的系数

$ \begin{array}{l} {\alpha _1} = {\xi _x}\frac{{{\partial ^2}x}}{{\partial {\xi ^2}}} + {\xi _y}\frac{{{\partial ^2}y}}{{\partial {\xi ^2}}},{\alpha _2} = 2\left( {{\xi _x}\frac{{{\partial ^2}x}}{{\partial \xi \partial \eta }} + {\xi _y}\frac{{{\partial ^2}y}}{{\partial \xi \partial \eta }}} \right),\\ {\alpha _3} = {\xi _x}\frac{{{\partial ^2}x}}{{\partial {\eta ^2}}} + {\xi _y}\frac{{{\partial ^2}y}}{{\partial {\eta ^2}}},{\alpha _4} = {\eta _x}\frac{{{\partial ^2}x}}{{\partial {\xi ^2}}} + {\eta _y}\frac{{{\partial ^2}y}}{{\partial {\xi ^2}}},\\ {\alpha _5} = 2\left( {{\eta _x}\frac{{{\partial ^2}x}}{{\partial \xi \partial \eta }} + {\eta _y}\frac{{{\partial ^2}y}}{{\partial \xi \partial \eta }}} \right),{\alpha _6} = {\eta _x}\frac{{{\partial ^2}x}}{{\partial {\eta ^2}}} + {\eta _y}\frac{{{\partial ^2}y}}{{\partial {\eta ^2}}}. \end{array} $ (3b)

由于采用静压分布,垂向动量方程忽略,垂向流速通过连续性方程求解。为了避免垂向的网格划分,这里垂向采用谱方法计算,假定水平流速分量的垂向分布能够由以相对高度为自变量的有限项正交多项式基函数近似表示,可以是Legendre多项式或Chebyshev多项式,这里采用Legendre多项式

$ \begin{array}{l} u \approx \sum\limits_{i = 0}^N {{c_i}{p_i}} ,v \approx \sum\limits_{i = 0}^N {{d_i}{p_i}} ,\\ {v_{\rm{z}}} \approx \sum\limits_{i = 0}^N {{e_{{\rm{z}}\;i}}{p_i}} ,{v_{\rm{h}}} \approx \sum\limits_{i = 0}^N {{e_{{\rm{h}}i}}{p_i}} . \end{array} $ (4)

其中:N为阶数,pii阶多项式,各项系数为待求变量。

1.2 加权余量法

为确定多项式系数,水平动量方程需要采用加权余量法方法进行变换,试函数与基函数相同。ξ方向动量方程转化为

$ \begin{array}{l} \int_{{z_{\rm{b}}}}^H {{u^{n + 1}}{p_k}{\rm{d}}z} = \int_{{z_{\rm{b}}}}^H {{u^{ * ;n}}{p_k}{\rm{d}}z} + \Delta t\int_{{z_{\rm{b}}}}^H {\left( {{f^n} + D_{\rm{h}}^n} \right){p_k}{\rm{d}}z} + \\ {p_k}\Delta t{v_{\rm{z}}}\frac{{\partial {u^n}}}{{\partial z}}\left| {_{{z_{\rm{b}}}}^H} \right. - \Delta t\int_{{z_{\rm{b}}}}^H {\frac{{\partial {p_k}}}{{\partial z}}{v_{\rm{z}}}\frac{{\partial {u^n}}}{{\partial z}}{\rm{d}}z} ,k = 0, \cdots ,N. \end{array} $ (5)

其中:上标n为前一时间步;*为经过对流项计算后的数值;f为压力项;Dh为水平紊动黏性项;这里积分变量独立于ξη。利用多项式的正交性质,各积分项均可化为多项式系数的函数,得到ξ方向N+1个关于多项式系数的待求解方程

$ \begin{array}{*{20}{c}} {\frac{h}{2}\sum\limits_{j = 1}^N {{c_j}^{n + 1}{B_{jk}}} = \frac{h}{2}\sum\limits_{j = 1}^N {{c_{j * }}^n{B_{jk}}} + \Delta t\frac{h}{2}{f^n}{B_{0k}} + }\\ {\frac{\partial }{{\partial \xi }}\left( {\sum\limits_{m = 0}^N {{e_{{\rm{h}}m}}\xi _r^2\frac{\partial }{{\partial \xi }}} \sum\limits_{j = 0}^N {{c_j}^n} } \right){C_{jmk}} + }\\ {\frac{\partial }{{\partial \eta }}\left( {\sum\limits_{m = 0}^N {{e_{{\rm{h}}m}}\eta _r^2\frac{\partial }{{\partial \eta }}} \sum\limits_{j = 0}^N {{c_j}^n} } \right){C_{jmk}} + \varphi \Delta t{\nu _{\rm{z}}}\frac{{\partial {u^n}}}{{\partial z}}\left| {_{{z_{\rm{b}}}}^H} \right. - }\\ {\Delta t\frac{2}{h}\sum\limits_{m = 0}^N {{e_{{\rm{z}}m}}} \sum\limits_{j = 0}^N {{c_j}^n{D_{jk,m}}} ,\;\;\;\;k = 0,1, \cdots ,N.} \end{array} $ (6)

其中与多项式有关的积分常数为

$ \begin{array}{l} {B_{jk}} = \int_{ - 1}^1 {{p_j}{p_k}{\rm{d}}\zeta } ,\\ {C_{jmk}} = \int_{ - 1}^1 {{p_k}{p_m}{p_j}{\rm{d}}\zeta } ,\\ {D_{jk,m}} = \int_{ - 1}^1 {\frac{{\partial {p_k}}}{{\partial \zeta }}{p_m}\frac{{\partial {p_j}}}{{\partial \zeta }}{\rm{d}}\zeta } . \end{array} $ (7)

水面剪切阻力假定为零。床面剪切力

$ {v_{\rm{z}}}\frac{{\partial u}}{{\partial z}}\left| {_{{z_{\rm{b}}}}} \right. = {U_ * }{u_ * } = {C_{\rm{d}}}{U_{\rm{b}}}{u_{\rm{b}}}. $ (8)

其中:Cd为床面剪切力系数;U*u*为总摩阻流速数值和ξ方向的协变分量;Ubub为近底床面流速数值和ξ方向的协变分量。加权余量法已将水面和床面的边界条件显式包含在求解方程中。

1.3 对流项

采用SLS求解具有双曲性质的对流项,已经得到深入研究[17-19]。SLS原理是对流项等价于单位时间内质点沿特征线方向逆向追踪前后的速度变化量。由于谱方法垂向不离散没有离散节点,计算对流项时需要根据谱方法配置点法(Collocation)的思想,由阶数确定垂向线上的Gauss积分点的垂向位置及流速,作为追踪过程的起点。由Gauss数值积分原理得

$ \begin{array}{l} \int_{ - 1}^1 {{p_i}J\left( {\frac{\partial }{{\partial \xi }}\frac{{uu}}{J} + \frac{\partial }{{\partial \eta }}\frac{{uv}}{J} + \frac{\partial }{{\partial z}}\frac{{uw}}{J}} \right){\rm{d}}\zeta } \approx \\ \sum\limits_{j = 1}^N {{{\left[ {{\omega _j}{p_i}J\left( {\frac{\partial }{{\partial \xi }}\frac{{uu}}{J} + \frac{\partial }{{\partial \eta }}\frac{{uv}}{J} + \frac{\partial }{{\partial z}}\frac{{uw}}{J}} \right)} \right]}_{\zeta = {\zeta _j}}}} . \end{array} $ (9)

其中:ωj为积分权系数;ζj为积分点。对流项积分转化为各个积分点对流项的加权和。在各积分点处采用SLS求解对流项

$ - J\left( {\frac{\partial }{{\partial \xi }}\frac{{uu}}{J} + \frac{\partial }{{\partial \eta }}\frac{{uv}}{J} + \frac{\partial }{{\partial z}}\frac{{uw}}{J}} \right) \approx \frac{{{u^{ * ;n}} - {u^n}}}{{\Delta t}}. $ (10)

其中:un为追踪起点处前一时间步(nth)计算的uu*; n为沿特征线逆向追踪Δt后当地u。逆向追踪过程小时间步长时可以采用单步追踪的立方插值拟质点(cubic-interpolated pseudo-particle,CIP)[18]方法, 大时间步长可采用多步追踪的Eulerian Lagrangian方法(ELM)[20]。相对于Eulerian方法显式求解对流项受到Courant数不能超过1的时间步长限制而言,SLS理论上自动迎风且不受Courant数的限制,Courant数仅仅作为时间步长的参考。

这里对流项求解采用多步法。选择Gauss积分点作为追踪起点的垂向相对位置。将时间步长Δt细分成m子步,u动量方程的追踪起点平面位置位于u的定义处,确定k子步时追踪位置Pk*,插值得到该点流速,确定k+1子步位置Pk+1*,重复该过程直到位置Pm*。最后加权积分得到动量方程的对流项。

1.4 水位求解

垂向从床面(z=zb)到水面(z=H)积分连续性方程(1),结合垂向流速的运动边界条件得到水位方程。

$ \frac{\partial }{{\partial t}}\left( {\frac{h}{J}} \right) + \frac{\partial }{{\partial \xi }}\left( {\frac{{h\bar u}}{J}} \right) + \frac{\partial }{{\partial \eta }}\left( {\frac{{h\bar v}}{J}} \right) = 0. $ (11)

其中uvuv的水深平均值。当水平流速分量求解得出后,水位通过式(11)求得。

1.5 数值算法

该模型建立在平面交错网格上,水位和流速协变分量的多项式系数为待求解变量。选用Legendre多项式,数值算法按以下顺序求解:

步骤1 每一时间步开始,显式计算水平动量方程中的扩散项和压力梯度项;

步骤2 半隐格式迭代计算水位和水平动量方程中的压力梯度项,直至水位达到给定的收敛标准;

步骤3 采用SLS计算对流项;

步骤4 采用连续性方程计算垂向流速。

2 急弯水槽试验验证

采用该模型对Blanckaert[3, 21]急弯水槽试验进行模拟。水槽平面形态如图 1所示,流速和水位测量断面依次为弯道上游M25(字母表示上下游,数字表示距离或圆心角)、M05,弯道段S15、S30、S60、S90、S120、S150、S180,弯道下游P05、P15、P25。弯道为平坡矩形水槽,宽B=1.3 m,弯道圆心角为193°,进口直线段9 m,出口直线段5 m,中轴线曲率半径R=1.7 m。选取恒定流工况水深0.159 m,流量0.089 m3/s的一组实验。侧壁为水力光滑的有机玻璃。模拟时不考虑侧壁阻力。计算区域进出口直线段均设置为2.6 m,出口水位0.151 5 m。离散网格数目纵向130 ×横向20。计算时间步长为0.01 s。设置多项式阶数取值为N=0~6,调整阻力系数使得进口水位与实测一致,然后分析流场信息。Zeng等[22]采用Reynolds时均(Reynolds-averaged Navier-Stokes, RANS)(记为Z-R)、van Balen等[23]采用大涡模拟(large eddy simulation, LES)(记为vB-L)和RANS模拟(记为vB-R)2种方法分别对该试验进行了模拟研究。

图 1 文[3]的试验水槽形态

水平和垂向涡黏性系数均采用抛物型分布,同时考虑紊动的各向异性,假定水平涡黏性系数是垂向的两倍, 则有

$ {\nu _{\rm{h}}} = 2{\nu _{\rm{z}}} = \frac{3}{2}{{\bar \nu }_{\rm{h}}}\left( {1.01 - {\zeta ^2}} \right), $ (12)
$ {{\bar \nu }_{\rm{h}}} = \frac{\kappa }{6}h{u_ * }. $ (13)

其中κ为Kármdn常数取值0.41。

不同工况下调整Cd数值,使得进口水位与实测一致,模拟运行200 s,水流均达到稳定,所需计算机单核CPU时间与多项式阶数的关系如表 1所示,其中N=0对应不考虑垂向变化的平面二维模型情况。从N=1开始,随着阶数增加,CPU时间呈线性增加。

表 1 水槽试验验证结果
阶数/其他模型 CPU用时/s 横向网格 垂向网格或Gauss点数 主流线位置均差/% S90中线流速均差/(m·s-1) S180中线流速均差/(m·s-1)
纵向 横向 纵向 横向
0 100 20 1 29.2
1 189 20 2 17.3
2 294 20 3 4.8 0.030 0.015 0.029 0.035
3 459 20 4 4.5 0.031 0.012 0.036 0.039
4 644 20 5 5.4 0.028 0.013 0.035 0.041
5 876 20 6 6.3 0.027 0.012 0.035 0.038
6 1 148 20 7 5.2 0.025 0.012 0.035 0.038
Z-R 101 35 14.7
vB-R 192 24 8.6
vB-L 192 24 6.6

图 2a为实测垂向平均流速,图 3a为实测水位分布。可以看到水流进入弯道后,水位出现横向比降来整体平衡弯曲流动离心力,在弯道进口处,内岸水位降低,外岸水位抬升,内岸水流加速形成主流线,靠近外岸侧出现低流速区。而二次流的横向动量输运使得沿程各断面主流线的位置逐渐向外岸偏移。在弯道出口,横向比降逐渐消失,内岸减速,外岸加速,水动力轴线贴近外岸。实测主流线位置从弯道进口断面的内岸沿程逐渐移动到出口断面的外岸。采用不同阶数的模型进行计算,模拟结果因阶数不同出现差异。图 2bN=0对应平面二维模型,没有二次流输运作用,弯道主流线沿程位置沿程不发生偏移,一直位于弯道内侧,仅在弯道出口处快速靠近外岸。图 2cN=1对应的流速垂向分布为线性,纵向流速受二次流影响仅调整速度垂向梯度和大小,模型不能充分考虑二次流对主流的反馈作用,主流线位置调整较快,在弯道前半部分即移动至外岸侧。N>1的模拟结果与实测值较为接近, 随着N增加变化不大,图 2d2e给出N=2和N=6的流速分布结果,均有效预测了主流线的调整过程。不同N值的水位模拟结果均和实测十分接近,图 3b给出了N=6的情况。

图 2 (网络版彩图)弯道垂向平均流速实测与模拟结果

图 3 (网络版彩图)水位平面分布图

选择水动力轴线的沿程位置比较与其他三维模型的模拟精度。图 4给出了弯道主流线沿程横向偏移的模拟与实测情况,纵轴表示各断面水动力轴线距离内岸的相对距离,内岸为0外岸为1,并将不同阶数N主流线的沿程横向摆动模拟结果与实测值的平均误差统计在表 1,其中Z-R由文[22]得到,vB-R由文[23]得到,两者的RANS模拟出的主流线位置在弯道内相对实测结果均偏向内岸,似乎黏着于内岸,而大涡模拟的主流线却十分贴近实测值。文[23]将主流线位置的差异模糊地归结于紊流模型和LES计算的不稳定性。尽管两者采用复杂的三维RANS和LES模型,并且在较为精细的网格上计算,模型在模拟主流线位置上的精度和本文的简化模型N>1时相当,N>1对应弯道内主流线相对位置的平均误差均在0.07以内,随着N值增加误差没有明显变化。通过本文给定的垂向涡黏性系数是水平涡黏性系数的一半,能够成功模拟出主流线位置的沿程变化,因此本文推断本质原因应该是紊流的各向异性。涡黏性系数的方向性差异是紊动强度在垂向和水平方向上的差异的表现。水平涡黏性系数比垂向涡黏性系数大,表示动量更容易在水平方向传递,弯道主流线也就容易发生摆动。为了论证这一观点,对涡黏性系数水平和垂向相等即各相同性的情况进行计算比较。设定多项式阶数为6,模拟的弯道垂向平均流速如图 2f所示,通过与考虑各向异性的图 2e对比,可以明显看出主流线位置沿程偏移量的差异。紊动各相同性计算的主流线位置相对于各向异性更靠近内岸,与图 4中采用各向同性紊流模型的RANS结果一致。因此可以确定在急弯水流条件下紊动各向异性控制了主流线的偏移过程。同时Blanckaert在文[3]中给出实测紊动能三方向分量差异支持了涡黏性系数的各向异性设定。

图 4 主流线横向相对位置的沿程变化

为了进一步评估N=2~6的模型模拟性能,需要对流速分布进行检验。图 5给出了断面S90和S180纵向流速等值线图,选择断面中心线作纵向流速和横向流速对比,如图 6图 7表 1统计了纵向和横向流速平均误差,模型模拟结果和实测值吻合良好。S90断面位于弯道中段,横向流速垂向平均值0.429 m/s,二次流已经完全发展。随着N值增加,纵向和横向流速的模拟精度略微改善,平均误差从0.03 m/s降低为0.025 m/s,精度一般。横向流速平均误差从0.015 m/s降低为0.012 m/s,精度均较好。S180断面靠近出口,横向流速垂向平均值0.448 m/s,流速分布的模拟精度受到上游误差的传播的影响而明显降低。纵向和横向流速的模拟结果较差,模拟横向流速在近床面比实测偏大,模拟纵向流速较实测偏小。总体来说模拟与实测的横向流速在分布的大小和形状上与实测值基本吻合,证明该模型具有三维模拟的能力。本模型中的垂向流速是由连续方程计算得来,垂向流速不满足动量守恒,近壁面区域模拟精度会明显下降。本模型采用网格较粗,在两个断面靠近内岸区域流动复杂,较大的流速位于底部,模拟的该区域流速偏小。该处区域采用非静压模型更为恰当。

图 5 (网络版彩图)断面纵向流速分布实测与N=6模拟对比

图 6 S90断面中心线纵向(右侧)与横向(左侧)流速垂向分布的模拟与实测值

图 7 S180断面中心线纵向(右侧)与横向(左侧)流速垂向分布的模拟与实测值

3 讨论

模型采用的多项式阶数有限,当N>2随着阶数增加,准确度没有明显提高,与谱方法的谱精度并不违背,误差出在模型采用的基本方程和模型假定上。当流速的垂向分布十分复杂时,低阶多项式已经不能满足精度要求,系统误差会很大。涡黏性是流速分布的另外一个控制因素,较大的流速梯度会带来较大涡黏性反过来抑制流速梯度。本算例中弯道下游段流速模拟不能很好反映实测结果,与涡黏性系数垂向分布粗略采用的抛物线型仅与当地河床阻力关联,而没有采用严密的封闭关系与流速进行耦合计算有关。同时测试算例中,模型没有考虑侧壁阻力,导致近壁面处的计算流速会偏大。

对于三维浅水方程,本文垂向不做网格离散,采用正交多项式拟合流速的垂向分布,通过加权余量法建立垂向谱方法并通过实例计算论证了这一求解途径是可行的。算例中的急弯水槽底部为平坡,较为理想,可以采用谱方法计算。对于实际河流的复杂边界条件,谱方法受到限制,需要采用适应性强的谱元法。通过采用SLS成功求解了动量方程的对流项,将流速的垂向分布计算转化成计算Legendre多项式积分零点处变化量的加权和,这种转化精度很高。在ζ < -1和ζ>1时并没有物理意义,通过水面的动量通量应该为零,但这在多项式阶数较低的情况下并不能严格满足,而是随着阶数增加呈现渐进性满足。床面阻力通过设置在床面上的非零流速计算,阻力系数不独立。急弯弯道靠近侧壁处采用静压假定误差较大,N=0作为一种特例,式(6)与η方向公式均可简化,联立式(11)和定解条件,构成垂向平均的二维明渠水动力学模型。

4 结论

本文提出的简化三维模型,采用平面二维与垂向谱方法结合,有很好的计算效率和精度。该模型通过对急弯明渠水流模拟,模拟的水流结构与试验结果匹配较好,并且低阶多项式的计算量与二维模型相当。模型采用零方程模型计算涡黏性系数并考虑各向异性,能够和复杂三维水动力学模型一样模拟出急弯的水流结构。在急弯水流流速的垂向分布、速度垂向平均值的平面分布、主流线位置上,都显示该模型能够进行很好的模拟急弯水流结构。

参考文献
[1]
ENGELUND F. Flow and bed topography in channel bends[J]. Journal of the Hydraulics Division, 1974, 100(3): 1631-1648.
[2]
KIKKAWA H, IKEDA S, KITAGAWA A. Flow and bed topography in curved open channels[J]. ASCE, 1976, 102(9): 1327-1342.
[3]
BLANCKAERT K. Saturation of curvature-induced secondary flow, energy losses, and turbulence in sharp open-channel bends:Laboratory experiments, analysis, and modeling[J]. Journal of Geophysical Research:Earth Surface, 2009, 114(F3): F03015.
[4]
BLANCKAERT K, DE VRIEND H J. Nonlinear modeling of mean flow redistribution in curved open channels[J]. Water Resources Research, 2003, 39(12): 1-6.
[5]
BLANCKAERT K, DE VRIEND H J. Meander dynamics:A nonlinear model without curvature restrictions for flow in open-channel bends[J]. Journal of Geophysical Research:Earth Surface, 2010, 115(F4): F04011.
[6]
HOSODA T, NAGATA N, KIMURA I, et al. A depth averaged model of open channel flows with lag between main flows and secondary currents in a generalized curvilinear coordinate system[J]. Advances in Fluid Modeling and Turbulence Measurements, 2002, 63-70.
[7]
ONDA S, HOSODA T, KIMURA I. Depth averaged model considering transformation of downstream velocity in curved channel in generalized curvilinear coordinate system and its verification[J]. Annual Journal of Hydraulic Engineering, 2006, 50: 769-774. DOI:10.2208/prohe.50.769
[8]
KIMURA I, ONDA S, HOSODA T, et al. Computations of suspended sediment transport in a shallow side-cavity using depth-averaged 2-D models with effects of secondary currents[J]. Journal of Hydro-environment Research, 2010, 4(2): 153-161. DOI:10.1016/j.jher.2010.04.008
[9]
YANG F, KIMURA I, SHIMIZU Y, et al. Computations of open channel flows in a shallow sharp bend using depth-averaged 2D models with secondary flow effects[J]. Journal of Japan Society of Civil Engineers, 2016, 72(2): 505-513.
[10]
BERNARD R S, SCHNEIDER M L. Depth-averaged numerical modeling for curved channels[R/OL].[2017-12-15]. https://www.researchgate.net/publication/235054674_Depth-Averaged_Numerical_Modeling_for_Curved_Channels.
[11]
UCHIDA T, FUKUOKA S. Numerical calculation for bed variation in compound-meandering channel using depth integrated model without assumption of shallow water flow[J]. Advances in Water Resources, 2014, 72: 45-56. DOI:10.1016/j.advwatres.2014.05.002
[12]
ODGAARD, J A. Meander flow model: I. Development[J]. Journal of Hydraulic Engineering, 112(12), 1117-1136, 1986.
[13]
IKEDA S, NISHIMURA T. Flow and bed profile in meandering sand-silt rivers[J]. Journal of Hydraulic Engineering, 1986, 112(7): 562-579. DOI:10.1061/(ASCE)0733-9429(1986)112:7(562)
[14]
JOHANNESSON H, Parker G. Secondary flow in mildly sinuous channel[J]. Journal of Hydraulic Engineering, 1989, 115(3): 289-308. DOI:10.1061/(ASCE)0733-9429(1989)115:3(289)
[15]
CANUTO C, HUSSAINI M Y, QUARTERONI A, et al. Spectral Methods in Fluid Dynamics[M]. New York: Springer Verlag, 1988.
[16]
KIMURA I, HOSODA T, ONDA S. Numerical simulator on full staggered boundary fitted coordinate system for the analysis of 3D turbulent flows in open channels[J]. Yokkaichi University Journal of Environmental and Information Sciences, 2002, 5(1-2): 145-170.
[17]
EWING R E, WANG H. A summary of numerical methods for time-dependent advection-dominated partial differential equations[J]. Journal of Computational and Applied Mathematics, 2001, 128(1): 423-445.
[18]
JANG C L, SHIMIZU Y. Numerical simulation of relatively wide, shallow channels with erodible banks[J]. Journal of Hydraulic Engineering, 2005, 131(7): 565-575. DOI:10.1061/(ASCE)0733-9429(2005)131:7(565)
[19]
YANG F, FU X. Eulerian-Lagrangian method based on streamline[C]//Proceeding of 2013 IAHR World Congress. Chengdu, China, 2013: 1-10.
[20]
ZHANG Y, BAPTISTA A M, MYERS E P. A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems:I. Formulation and skill assessment[J]. Continental Shelf Research, 2004, 24(18): 2187-2214. DOI:10.1016/j.csr.2004.07.021
[21]
BLANCKAERT K. Flow and turbulence in sharp open-channel bends[D]. Switzerland: Ecole Polytechnique Fe'de'rale Lausanne, 2002.
[22]
ZENG J, CONSTANTINESCU G, BLANCKAERT K, et al. Flow and bathymetry in sharp open-channel bends:Experiments and predictions[J]. Water Resource Research, 2008, 44(9).
[23]
VAN BALEN W, BLANCKAERT K, UIJTTEWAAL W S J. Analysis of the role of turbulence in curved open-channel flow at different water depths by means of experiments, LES and RANS[J]. Journal of Turbulence, 2010(11): N12.