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

一类二次流修正方法^[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为重力加速度；x、y和z为直角坐标；ξ、η为平面曲坐标；ν_z、ν_h为垂向与水平涡黏性系数；J=ξ_xη_y-ξ_yη_x为Jacobian行列式；u、v为水平流速在曲坐标系下的协变分量；w为垂向流速；ξ_r、η_r为曲坐标系两个方向Lame系数的倒数；ζ=2(z-z_b)/h-1为相对高度；z_b为床面高程。动量方程中的水平扩散项表达式已做近似简化处理，压力梯度项系数

$ \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为阶数，p_i为i阶多项式，各项系数为待求变量。

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为压力项；D_h为水平紊动黏性项；这里积分变量独立于ξ和η。利用多项式的正交性质，各积分项均可化为多项式系数的函数，得到ξ方向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)

其中：C_d为床面剪切力系数；U_*、u_*为总摩阻流速数值和ξ方向的协变分量；U_b、u_b为近底床面流速数值和ξ方向的协变分量。加权余量法已将水面和床面的边界条件显式包含在求解方程中。

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)

其中：uⁿ为追踪起点处前一时间步(n^th)计算的u；u^{*; 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子步时追踪位置P_k^*，插值得到该点流速，确定k+1子步位置P_k+1^*，重复该过程直到位置P_m^*。最后加权积分得到动量方程的对流项。

1.4 水位求解

垂向从床面(z=z_b)到水面(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)

其中u、v为u、v的水深平均值。当水平流速分量求解得出后，水位通过式(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 m³/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种方法分别对该试验进行了模拟研究。

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

$ {\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。

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

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

选择水动力轴线的沿程位置比较与其他三维模型的模拟精度。图 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]中给出实测紊动能三方向分量差异支持了涡黏性系数的各向异性设定。

为了进一步评估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，流速分布的模拟精度受到上游误差的传播的影响而明显降低。纵向和横向流速的模拟结果较差，模拟横向流速在近床面比实测偏大，模拟纵向流速较实测偏小。总体来说模拟与实测的横向流速在分布的大小和形状上与实测值基本吻合，证明该模型具有三维模拟的能力。本模型中的垂向流速是由连续方程计算得来，垂向流速不满足动量守恒，近壁面区域模拟精度会明显下降。本模型采用网格较粗，在两个断面靠近内岸区域流动复杂，较大的流速位于底部，模拟的该区域流速偏小。该处区域采用非静压模型更为恰当。

3 讨论

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

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

4 结论

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