浅水波方程混合旋转通量旋转角的研究

Study on the Rotation Angle of Mixed Rotational Flux in the Shallow Water Equation

摘要: 本文在浅水波方程的混合旋转通量的基础上，对通量旋转角进行分析及研究。在单元界面法向方向采用具有良好鲁棒性且可消除红斑现象的HLL格式，在单元界面切线方向上采用满足热力学第二定律的熵稳定格式，两者进行结合，时间方向采用三阶强稳定Runge-Kutta方法，得到混合旋转通量算法。基于Euler方程通过压力加权函数定义旋转角度的策略，对浅水波方程混合旋转通量的旋转角给出不同的方案并进行数值实验，验证算法的有效性，确定最佳旋转角度。

Abstract: Based on the mixed rotational flux of the shallow water equations, this paper conducts an analysis and research on the flux rotation angle. For the normal direction of the element interface, the HLL scheme—characterized by excellent robustness and the ability to eliminate the “red spot” phenomenon is adopted; for the tangential direction of the element interface, an entropy-stable scheme that satisfies the second law of thermodynamics is used. By combining these two schemes and applying the third-order strongly stable Runge-Kutta method in the time direction, a mixed rotational flux algorithm is derived. Drawing on the strategy of defining the rotation angle via a pressure-weighted function based on the Euler equations, different schemes for the rotation angle of the mixed rotational flux in the shallow water equations are proposed. Numerical experiments are carried out to verify the effectiveness of the algorithm and determine the optimal rotation angle.

1. 引言

计算流体力学是借助数值方法与计算机，求解描述流体运动的控制方程以模拟、分析流体流动现象的学科，它通过将连续流场离散为网格节点，把偏微分方程转化为代数方程组迭代求解，从而获取速度、压力等物理量分布，广泛用于航空航天、汽车设计等领域以降低实验成本、研究复杂流动。计算流体力学的核心任务之一就是研发稳定、高精度的数值方法，求解这类方程并准确捕捉激波、间断等复杂流动结构，确保数值解既符合数学稳定性要求，又能反映真实的流体物理规律。浅水波方程、Euler方程等均可以写成双曲守恒律方程[1]的形式。双曲守恒律方程是计算流体力学的重要数学基础，目前求解双曲守恒律方程的数值方法研究围绕精度提升、间断捕捉能力强化、计算效率优化及复杂场景适应性展开，形成了传统基础方法与高阶高精度方法并行发展的格局。早期发展的经典数值方法有有限差分法[2] [3]及有限体积法，为后续方法奠定了理论基础，有效减少了数值耗散与波动；TVD格式[4]作为高分辨率格式的重要分支，实现了间断附近的无振荡求解，能在保持一阶精度捕捉间断的同时，在光滑区域达到高阶精度，鲁棒性强，有效平衡了数值稳定性与间断分辨率；本质无振荡(ENO)方法[5]与加权本质无振荡(WENO)方法[6] [7]是高阶高精度方法的主流方向，ENO方法解决了传统方法在间断附近的振荡问题，WENO方法在此基础上引入非线性权重，进一步提升了精度与稳定性。国内外学者通过优化权重函数、改进插值多项式构造，发展出高阶HWENO格式等变体，高阶精度与高分辨率兼具，能同时处理光滑区域与强间断。

对双曲守恒律方程通量函数的离散方向，Levy等[8]研究了多维控制方程的旋转黎曼求解器，为了更好地捕捉激波，同时保留耗散通量的鲁棒性，研究者采用混合格式。Hiroaki Nishikawa与Keiichi Kitamura [9]提出一种基于旋转Rieman求解器的混合格式，核心是将Roe格式与HLL格式耦合，通过自适应选择求解方向，在激波法向自动应用少波求解器抑制红斑现象，在剪切层区域应用全波求解器避免过度耗散。封建湖、刘友琼、任炯团队[10]将HLL格式与旋转通量思想结合，创新地在激波法向采用HLL格式抑制“红斑现象”，在激波传播方向采用HLLC格式减少过度耗散，通过旋转Riemann求解器将二维问题转化为类一维处理，解决了传统HLL类格式在多维强激波问题中易出现的非物理振荡，格式结构简洁且边界层分辨率优异。贾豆，郑素佩等[11]将熵稳定格式与HLL格式结合提出了一种旋转混合格式，用于求解二维欧拉方程，此格式具有数值稳定、分辨率高等优点。

本文基于混合旋转通量的基础上，对通量旋转角进行分析及研究。介绍了二维浅水波方程、HLL-熵稳定混合旋转格式[12]等基础知识，将压力加权思想运用到求解浅水波方程的过程中，对浅水波方程混合旋转通量的旋转角给出不同的量值进行尝试，寻找性能良好的旋转角。

2. 控制方程

二维浅水波方程[12]为

$U_t+F{\left(U\right)}_x+G{\left(U\right)}_y=-S,$ (1)

其中，

$U=\left[\begin{array}{c}h\\ hu\\ hv\end{array}\right],F\left(U\right)=\left[\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\\ huv\end{array}\right],G\left(U\right)=\left[\begin{array}{c}hv\\ huv\\ h{v}^2+\frac{1}{2}g{h}^2\end{array}\right],$

$h$ 为水深，$g$ 为重力加速度，x方向上的水速$u$ ，y方向上的水速$v$ ，源项$S$ 。其中$S$ 可分为摩擦项和倾斜项，

$S=\left[\begin{array}{c}0\\ gh\frac{\partial z_b}{\partial x}+gh{S}_{f_x}\\ gh\frac{\partial z_b}{\partial y}+gh{S}_{f_y}\end{array}\right],$

$z_b$ 表示底部地势函数，$gh\left(\partial z_b/\partial x\right)$ 和$gh\left(\partial z_b/\partial y\right)$ 表示水下底部作用力沿x方向和y方向上的分力，$gh{S}_{f_x}$ 和$gh{S}_{f_y}$ 为水下底面摩擦力沿x方向和y方向的量，$S_{f_x}$ 和$S_{f_y}$ 为x方向和y方向上的摩阻比率，有

$S_{f_x}=\frac{u}{K^2}\frac{{\left(u^2+v^2\right)}^{1/2}}{h^{4/3}},S_{f_y}=\frac{v}{K^2}\frac{{\left(u^2+v^2\right)}^{1/2}}{h^{4/3}},$

$K$ 表示摩擦因数，通常情况取$K=50$ 。

3. HLL-熵稳定混合旋转格式[13]

为了使数值结果具有更好的鲁棒性和更高的分辨率，本文基于浅水波方程通量函数的旋转不变性，结合熵稳定通量和HLL通量，构造一种旋转混合通量。根据浅水波方程的旋转不变性，将界面的单位外法线方向$n$ 分解，混合格式将$n$ 分解到两个正交方向$n_1,n_2$ 上，也就是选择好方向$n_1$ (本文选取的方向与x方向相同)，方向$n_2$ 必须是垂直于方向$n_1$ 的，即

$n_1\cdot n_2=0$ ，其中$\left|n_1\right|=\left|n_2\right|=1$

单元界面外法向量$n$ 投射到两个正交方向$n_1,n_2$ 上有

$n={\alpha }_1{n}_1+{\alpha }_2{n}_2,$

其中“$\cdot$ ”表示向量的点积，${\alpha }_1\ge 0,{\alpha }_2\ge 0$ ，且${\alpha }_1+{\alpha }_2=1$ 。

在$n_1,n_2$ 上构造熵稳定通量和HLL通量耦合成的混合通量格式，得

$H^{混}={\alpha }_1{H}^{HLL}+{\alpha }_2{H}^{ES}$ (2)

$H^{HLL}=\frac{S_R^{+}{H}_k\left(U_L\right)-S_L^{-}{H}_k\left(U_R\right)}{S_R^{+}-S_L^{-}}+\frac{S_R^{+}{S}_L^{-}}{S_R^{+}-S_L^{-}}\Delta U,$ (3)

$H_k=\left(F,G\right)\cdot n\left(k=x,y\right),n=\left(n_x,n_y\right),\Delta U=U_R-U_L,$

$S_R^{+}=\max \left(0,S_R\right),S_L^{-}=\min \left(0,S_L\right),$

$H^{ES}=H^{EC}-\frac{1}{2}R\Lambda R^{\text{T}}dV.$ (4)

其中$H^{EC}$ 为熵守恒通量，$\Lambda$ 为对角矩阵，$R$ 为右特征向量矩阵。利用算术平均进计算$H^{EC}$ 中的通量，则

$H^{EC}=n_{x}F+n_{y}G,$ (5)

$F=\left[\begin{array}{c}\frac{h_L+h_R}{2}\frac{u_L+u_R}{2}\\ \frac{h_L+h_R}{2}{\left(\frac{u_L+u_R}{2}\right)}^2+\frac{g}{4}\left(h_L^2+h_R^2\right)\\ \frac{h_L+h_R}{2}\frac{u_L+u_R}{2}\frac{v_L+v_R}{2}\end{array}\right],G=\left[\begin{array}{c}\frac{h_L+h_R}{2}\frac{v_L+v_R}{2}\\ \frac{h_L+h_R}{2}\frac{u_L+u_R}{2}\frac{v_L+v_R}{2}\\ \frac{h_L+h_R}{2}{\left(\frac{v_L+v_R}{2}\right)}^2+\frac{g}{4}\left(h_L^2+h_R^2\right)\end{array}\right].$

$H^{HLL}$ 为(3)式HLL数值通量，$H^{ES}$ 为(4)式熵稳定通量。

4. 旋转角的压力加权旋转策略[14]

针对Euler方程，密度项在数值通量中的额外耗散将抑制激波的不稳定性。同时，在粘性流动的数值模拟中，格式的精度也有所下降。由于耗散项含有旋转角且耗散具有单调性，进而找到一种改进旋转迎风格式性能的策略。这个简单的策略就是通过压力来控制旋转角度。

该策略的关键是通过局部激波强度来调节耗散，因此需要一种激波检测方法。Quirk提出了一种基于压力的转换函数用于HLLE格式的自适应应用。在这里，Kim等人的压力加权函数被用作激波强度的基本指标。其表达式为

$f=1-\min {\left(\frac{p_L}{p_R},\frac{p_R}{p_L}\right)}^3,$ (6)

$p_L,p_R$ 分别表示单元界面左右两边的压力，该函数值随界面两侧压差的增大而逐渐变化，且该函数是稳定的。

当需要计算通过某个单元界面的数值通量时，激波检测函数应检查所有相邻的界面。这里定义了一个函数，用于选择单元i和j交界面上的最大函数值$f_k$ ：

$\omega =\underset{k}{\max }\left(f_k\right),$ (7)

在光滑流动区域，各界面的$f_k$ 值大小相近，因此函数最大值对计算结果影响不大。

为简单起见，在一阶和二阶空间精度应用程序中函数$f$ 中的左右压力一般直接使用左右单元格中心值。在强激波区域，左右压差较大，此时函数值$\omega$ 约等于1，而在光滑流动区域，函数值会很小。因此，函数$\omega$ 适用于识别激波的存在。

将旋转角度定义为

$\alpha =\frac{\pi }{4}\omega .$ (8)

旋转角度取决于局部压差的大小。在强激波附近，旋转角约等于45度。相反，在没有强压力且不连续性的情况下，可以采用网格对齐或“近似”网格对齐的方案，此时旋转角度近似等于0度。

针对二维浅水波方程，方程中不存在压力的变量，基于压力加权方法中旋转角度的定义式，尝试在求解浅水波方程时旋转的最佳角度，故分别取旋转角$\alpha$ 为角$\frac{\pi }{\text{3}},\frac{\pi }{4},\frac{\pi }{\text{6}}$ ，运用到二维浅水波方程的旋转通量中，通过数值实验对比验证。

5. 数值算例

本节采用所构造的格式求解数值算例，取网格数均为$\text{10}0\times \text{10}0$ ，并对所得结果进行分析、比较。

例1 偏心圆形水柱演化

此问题的计算区域为$\left[-\text{10},\text{4}0\text{m}\right]\times \left[10,60\text{m}\right]$ ，当$S=0$ 时，方程(1)满足如下条件，

$\begin{array}{l}h\left(x,y,0\right)=\{\begin{array}{l}\text{8},\sqrt{{\left(x-\text{1}5\right)}^2+{\left(y-\text{3}5\right)}^2}<\text{8}\\ 1,\text{others}\end{array},\\ u\left(x,y,0\right)=v\left(x,y,0\right)=0,z_b\left(x,y\right)=0.\end{array}$

其中$CFL=0.\text{8}$ ，时间$t=\text{5},g=0.98\text{1}$ ，边界条件采用周期性条件。图1为偏心圆形水柱水深等值线图；图1(a)为偏心圆形水柱三维立体模拟图，图1(b)显示的是旋转通量中旋转角$\alpha =\pi /4$ 密度结果图，图1(c)显示的是旋转通量中旋转角$\alpha =\pi /3$ 密度结果图。图1(d)显示的是旋转通量中旋转角$\alpha =\pi /6$ 密度结果图。通过结果图对比可以发现旋转角$\alpha =\pi /4$ 时的结果对称性均良好，结构清晰，具有更高的分辨率。

(a) (b)

(c) (d)

Figure 1. Contour map of water depth for eccentric circular water column. (a) Three-dimensional stereogram; (b) Density of rotation angle $\alpha =\pi /4$ ; (c) Density of rotation angle $\alpha =\pi /3$ ; (d) Density of rotation angle $\alpha =\pi /6$

图1. 偏心圆形水柱水深等值线图。(a) 三维立体图；(b) 旋转角$\alpha =\pi /4$ 密度图；(c) 旋转角$\alpha =\pi /3$ 密度图；(d) 旋转角$\alpha =\pi /6$ 密度图

例2 干河床上的圆形冲击

此问题计算区域为$\left[-\text{15},\text{35}\text{m}\right]\times \left[\text{0},\text{50}\text{m}\right]$ ，当$S=0$ 时，方程(1)满足如下条件，

$\begin{array}{l}h\left(x,y,0\right)=\{\begin{array}{l}\text{15},\sqrt{{\left(x-\text{10}\right)}^2+{\left(y-\text{2}5\right)}^2}<\text{5}\\ 1,\text{others}\end{array},\\ u\left(x,y,0\right)=v\left(x,y,0\right)=0,z_b\left(x,y\right)=0.\end{array}$

其中$CFL=0.\text{8}$ ，时间$t=\text{3},g=0.98\text{1}$ ，边界条件采用周期性条件。图2为干河床上的圆形冲击问题水深等值线图。图2(a)为干河床上的圆形冲击三维立体模拟图，图2(b)显示的是旋转通量中旋转角$\alpha =\pi /4$ 密度结果图，图2(c)显示的是旋转通量中旋转角$\alpha =\pi /3$ 密度结果图，图2(d)显示的是旋转通量中旋转角$\alpha =\pi /6$ 密度结果图。通过图形对比可以看出，相较于其他角度，旋转角$\alpha =\pi /4$ 时对称性更好，分辨率更高，数值结果没有出现振荡或过渡抹平现象。

(a) (b)

(c) (d)

Figure 2. Contour map of water depth for circular impact on a dry riverbed. (a) Three-dimensional stereogram; (b) Density of rotation angle $\alpha =\pi /4$ ; (c) Density of rotation angle $\alpha =\pi /3$ ; (d) Density of rotation angle $\alpha =\pi /6$

图2. 干河床上的圆形冲击问题水深等值线图。(a) 三维立体图；(b) 旋转角$\alpha =\pi /4$ 密度图；(c) 旋转角$\alpha =\pi /3$ 密度图；(d) 旋转角$\alpha =\pi /6$ 密度图

例3 环形水柱演化问题

此问题计算区域为$\left[\text{0},\text{80}\text{m}\right]\times \left[\text{0},\text{80}\text{m}\right]$ ，当$S=0$ 时，方程(1)满足如下条件

$\begin{array}{l}h\left(x,y,0\right)=\{\begin{array}{l}\text{10},\text{12}\le \sqrt{{\left(x-\text{40}\right)}^2+{\left(y-\text{40}\right)}^2}<\text{20}\\ 1,\text{others}\end{array},\\ u\left(x,y,0\right)=v\left(x,y,0\right)=0,z_b\left(x,y\right)=0.\end{array}$

取$CFL=0.\text{8}$ ，时间$t=\text{4},g=0.98\text{1}$ ，边界条件采用周期性条件。图3表示的是环形水柱演化问题水深等值线图。图3(a)为环形水柱演化三维立体模拟图，图3(b)显示的是旋转通量中旋转角$\alpha =\pi /4$ 密度结果图，图3(c)显示的是旋转通量中旋转角$\alpha =\pi /3$ 密度结果图，图3(d)显示的是旋转通量中旋转角$\alpha =\pi /6$ 密度结果图。通过图形对比可以看出旋转角$\alpha =\pi /4$ 时得到的数值结果具有更高的分辨率，对称性好，没有出现非物理振荡。

(a) (b)

(c) (d)

Figure 3. Contour map of water depth for the evolution of an annular water column. (a) Three-dimensional stereogram; (b) Density of rotation angle $\alpha =\pi /4$ ; (c) Density of rotation angle $\alpha =\pi /3$ ; (d) Density of rotation angle $\alpha =\pi /6$

图3. 环形水柱演化问题水深等值线图。(a) 三维立体图；(b) 旋转角$\alpha =\pi /4$ 密度图；(c) 旋转角$\alpha =\pi /3$ 密度图；(d) 旋转角$\alpha =\pi /6$ 密度图

例4 二维圆形浅水波传播问题

此问题计算区域为$\left[-\text{3},\text{3}\text{m}\right]\times \left[-\text{3},\text{3}\text{m}\right]$ ，底部地势函数和初始条件为

$\begin{array}{l}{z}_b\left(x,y\right)=1+0.0\text{2}\cos \left(\pi x/\text{3}\right)\cos \left(\pi y/\text{3}\right)\\ h\left(x,y,0\right)=z_b\left(x,y,0\right)+\zeta \left(x,y,0\right),u\left(x,y,0\right)=v\left(x,y,0\right)=0\end{array}$

其中初始水深表达式为：$\zeta \left(x,y,0\right)=0.1\text{5}\exp \left(-\text{20}{R}^2\right),R=\sqrt{x^2+y^2}$ 。

取$CFL=0.\text{45}$ ，时间$t=\text{0}.\text{4},g=0.98$ ，边界条件采用周期性条件。图4表示的是二维圆形浅水波传播问题水深等值线图。图4(a)为二维圆形浅水波传播三维立体模拟图，图4(b)显示的是旋转通量中旋转角$\alpha =\pi /4$ 密度结果图，图4(c)显示的是旋转通量中旋转角$\alpha =\pi /3$ 密度结果图，图4(d)显示的是旋转通量中旋转角$\alpha =\pi /6$ 密度结果图。通过图形对比可以看出旋转角$\alpha =\pi /4$ 时得到的结果对称性均良好，结构清晰，没有出现数值振荡或过渡抹平的现象。

(a) (b)

(c) (d)

Figure 4. Water depth contours for 2D circular shallow water wave propagation. (a) Three-dimensional stereogram; (b) Density of rotation angle $\alpha =\pi /4$ ; (c) Density of rotation angle $\alpha =\pi /3$ ; (d) Density of rotation angle $\alpha =\pi /6$

图4. 二维圆形浅水波传播问题水深等值线图。(a) 三维立体图；(b) 旋转角$\alpha =\pi /4$ 密度图；(c) 旋转角$\alpha =\pi /3$ 密度图；(d) 旋转角$\alpha =\pi /6$ 密度图

6. 总结

本文在旋转通量混合格式基础上受压力加权方法的启发，对通量函数中的旋转角进行了进一步研究，结果表明二维浅水波方程无论在源项为零，还是非零源项的情况下，数值结果均在旋转角$\alpha =\frac{\pi }{4}$ 时结构更清晰，对称性更加良好，没有出现伪振荡或过渡抹平现象。

基金项目

2023年度宁夏师范学院校级科研重点D类项目(XJZDD2302)。