[发明专利]一种求解浅水问题模拟间断水流数值的方法

摘要

本发明公开了一种求解浅水问题模拟间断水流数值的方法，涉及计算流体力学领域，具体公开了一种求解二维双曲型方程组黎曼问题的数值方法，将其应用于求解浅水波方程，进行间断水流的数值模拟。本发明在二维无结构三角形网格中，求解双曲型浅水方程的黎曼问题，提高了数值精度和分辨率，探讨了间断水流的水动力特征。本发明基于无结构三角形网格建立了能够适应复杂地形和复杂计算区域中带有大梯度或间断水流现象的高分辨率数值模型，该方法可以较大程度的减小误差，提高模拟水动力特征的分辨率。

主权项

一种求解浅水问题模拟间断水流数值的方法，其特征在于，建立模拟间断水流问题的高分辨率数值模型，利用所述模型模拟地形和计算域中带有大梯度解或间断的水流现象，具体步骤如下：步骤一、建立计算网格，将物理计算区域进行无结构三角形网格分割，得到复数个三角形网格单元，将单一的网格单元作为控制元，其中第i个三角形控制元记为Ω_i；步骤二、通过守恒型物理量反映每个计算网格的水动力特征，在笛卡尔直角坐标系下，将平面二维浅水运动守恒形式的方程表示为：U_t+F(U)_x+G(U)_y＝S    (1)；其中，U为守恒量向量；F，G为通量向量；S为源项向量，t为时间，x和y为笛卡尔坐标系坐标，U_t,F(U)_x,G(U)_y分别表示为时间t和空间x,y方向的偏导数，具体表示为：U＝[h,uh,vh]^T；F＝[uh,u²h+gh²/2,uvh]^T；G＝[vh,uvh,v²h+gh²/2]^T；S＝[0,‑ghb_x,‑ghb_y]^T；其中，T表示向量的转置；u、v分别是x、y方向沿水深平均的流速分量,并且为关于x,y,t的函数；h是水深；g是重力加速度；b_x、b_y分别为河床比降在x、y方向上的分量；求解上述方程的黎曼解，对公式(1)，令向量函数E＝(F,G)，为梯度算子，将式(1)表示为：U_t+▽·E＝S    (2)；采用步骤一中划分的三角形网格单元对计算区域进行离散，在Ω_i上对方程组(2)进行积分，并按逆时针利用Green公式将面积分化为线积分，得：$\frac{d}{\mathrm{dt}}{\∫}_{{\mathrm{\Ω}}_i}\mathrm{Ud\Ω}=-\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{\∫}_{L_m}E\·n_m\mathrm{dl}+{\∫}_{{\mathrm{\Ω}}_i}\mathrm{Sd\Ω}---\left(3\right);$其中，为时间方向的常微分算子，dΩ是面积分微元，dl是线积分微元，L_m表示第i个三角形单元Ω_i的第m条边，m＝1,2,3,n_m＝(n_mx,n_my)＝(cosθ,sinθ)为L_m的单位外法向量，n_mx,n_my为n_m在x、y方向的分量,θ为外法线方向与水平方向的夹角；定义单元平均值：${\stackrel{\‾}{U}}_i=\frac{1}{\left|{\mathrm{\Ω}}_i\right|}{\∫}_{{\mathrm{\Ω}}_i}\mathrm{Ud\Ω},{\stackrel{\‾}{S}}_i=\frac{1}{\left|{\mathrm{\Ω}}_i\right|}{\∫}_{{\mathrm{\Ω}}_i}\mathrm{Sd\Ω},$其中，|Ω_i|为网格单元面积；利用两点Gauss求积公式求解式(3)中的线积分：其中，|L_m|为三角形边界线L_m的长度，G_q为三角形单元边L_m上的Gauss点，为相应的权系数，q为高斯点编号，构造一个数值通量格式近似替代E(U(G_q,t))，将三角形网格单元边界上内部和外部的变量值分别记为U_L、U_R，将通过Gauss点G_q的通量记为其中，U_L(G_q,t)，U_R(G_q,t)分别为U在Gauss点三角形网格单元内部和外部的值，由此得出控制方程(3)的有限体积半离散化格式：$\frac{d}{\mathrm{dt}}{\stackrel{\‾}{U}}_i=-\frac{1}{\left|{\mathrm{\Ω}}_i\right|}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}\left|L_m\right|\underset{q=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\ω}}_q\hat{E}\left(U_R(G_q,t)U_L(G_q,t)\right)\·n_m+{\stackrel{\‾}{S}}_i---\left(5\right);$步骤三、求解三角形网格单元边界上的HLL数值通量步骤四、采用三阶Runge‑Kutta时间离散方法对时间方向上的离散进行处理，得出三阶精度的时间离散格式；步骤五、采用三阶WENO格式重构方法求解U_L(G_q,t)，U_R(G_q,t)；步骤六、对每个网格单元，利用当前时间层上每个单元的守恒量向量积分平均值，求解下一时刻的守恒量向量积分平均值：经过步骤五对流体方程时间方向上的的离散后，再对流体方程进行空间离散，表达式为：$\frac{\mathrm{dU}}{\mathrm{dt}}=L\left(U\right);$其中，L(U)是空间离散化算子；设定第n个时间层上的U值为Uⁿ，则经过一个时间步长Δt后，下个时间层的近似值Uⁿ⁺¹表示为：U⁽¹⁾＝Uⁿ+ΔtL(Uⁿ)；$U^{\left(2\right)}=\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\left(1\right)}+\frac{1}{4}\mathrm{\ΔtL}\left(U^{\left(1\right)}\right);$$U^{n+1}=\frac{1}{3}{U}^n+\frac{2}{3}{U}^{\left(2\right)}+\frac{2}{3}\mathrm{\ΔtL}\left(U^{\left(2\right)}\right);$其中，U⁽¹⁾、U⁽²⁾为中间过渡向量。