下载:
-
水力系统中部分位置为明流, 部分位置为满流, 即明满流共存；或者同一位置不同时刻可能出现明流也可能出现满流, 即明满流交替出现, 这些都属于明满流过渡问题。水工隧洞、城市排水管道、水电站尾水系统中可能出现明满流过渡问题。另外, 严寒地区河道冰塞以及固定冰盖下电站调峰时亦可出现明满流过渡问题。明满流过渡会导致水力系统压力瞬间发生明显变化, 对系统的安全运行构成较大威胁。对明满流过渡过程伴随的压力变化进行准确预报是水力系统设计、运行管理的关键性技术问题。然而在数值模拟中, 明流和有压流的控制方程不同, 明满流过渡模拟存在较大困难。明满流过渡模拟方法主要有：激波跟踪法、刚性水体法和激波捕捉法。其中, 以Preissmann窄缝法(Preissmann Slot Method, PSM)为代表的激波捕捉法因其方法简单而广受欢迎, 著名的SWMM(Storm Water Management Model)模型引入了该方法[1]。PSM是在密闭管道顶部加窄缝, 明流和有压流就可以采用相同的控制方程, 无需人为制定满流的判别标准和记录流动状态的转换, 编程更为简单。但该方法的缺点是无法模拟负压和波速突变产生明显非物理振荡。
模拟负压方面, Vasconcelos等[2]提出了TPA模型, Kerger等[3]提出了负Preissmann窄缝(Negative Preissmann Slot)概念及方法, 二者都是根据通气条件来修改动量方程中的压力项, 使其具备负压模拟功能, 但会引起非物理振荡。在抑制非物理振荡方面主要有构建高分辨率格式和增加数值黏性两大类方法。高分辨率格式构建方面, 如TVD(Total Variation Diminishing)格式、TVB(Total Variation Bounded)格式、ENO(Essentially non-Oscillatory)格式、PPM(Piecewise Parabolic Method)格式等。这些格式本身带有一定的数值黏性, 计算过程中存在耗散, 对线性间断的模拟效果不好[4]。增加数值黏性抑制非物理振荡主要有3种方法：①滤波法, 类似滑动平均, 有效抑制非物理振荡, 但也可能将实际的物理振荡平滑掉, 计算精度不高；②数值扩散法, 类似Lax-Friedrichs方法, 引入数值扩散可有效抑制非物理振荡, 但使自由表面流动具备较小的Courant数, 影响计算精度；③波速转换法, 将自由表面流动波速平稳转换至压力波波速, 避免产生非物理振荡。波速转换法应用相对广泛, 例如采用变缝宽方法来抑制非物理振荡[5]；基于HLL近似Riemann求解器, 采用修正的阶梯公式计算波速, 从而抑制振荡[5-7]。将波速转换法与数值扩散法进行组合, 称混合通量法(Hybrid Flux Approach)。如将Roe近似Riemann求解器中的特征值表达式进行修正, 改变强间断所在单元的扩散性来抑制振荡[8]。将HLL近似Riemann求解器与中心通量求解器(FORCE)联合使用, 当水流接近管顶时, 保证水流具有较强的扩散性, 达到抑制非物理振荡的目的[9]。
采用PSM法模拟明满流过渡过程, 核心问题是抑制非物理振荡。在以往的研究中, 均通过增加数值黏性来抑制振荡, 人为因素大, 改变了水流运动的物理特性。另外, 近似Riemann求解器存在诸多问题, Roe线性化Riemann求解器在跨临界流模拟中必须进行熵修正, 避免出现非物理解；HLL简化Riemann求解器对接触间断模拟精度较差。精确Riemann求解器是非线性、全波族的求解器, 计算耗时仅比近似Riemann求解器高25%[10]。本文基于Godunov格式, 采用精确Riemann求解器, 在非物理振荡成因分析的基础上, 提出基于精确Riemann解的变量空间重构方法, 表达激波间断在单元内的空间分布状态, 以期从机理上抑制非物理振荡。
-
一维明渠非恒定流动可用双曲型偏微分方程来描述, 通常称为圣维南方程。守恒型式的圣维南方程组如下[11]：
$$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} = \mathit{\boldsymbol{S}}\\ \mathit{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{l}} A\\ Q \end{array}} \right]\quad \mathit{\boldsymbol{F}} = \left[ {\begin{array}{*{20}{c}} Q\\ {{Q^2}/A + g{I_1}} \end{array}} \right]\quad \mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{c}} 0\\ {gA({S_0} - {S_{\rm{f}}}) + g{I_2}} \end{array}} \right] \end{array} $$ (1) 式中: U为变量；F为通量；S为源项；t为时间；x为空间坐标；A=A(x, t)为过流断面面积；Q为流量；g为重力加速度；S0为床面比降；Sf为摩阻比降；I1和I2分别为静力矩和侧压力, $I_{1}=\int_{0}^{h}(h-\eta) b(x, \eta) \mathrm{d} \eta$, $I_{2}=\int_{0}^{h}(h-\eta) \frac{\partial b(x, \eta)}{\partial x} \mathrm{~d} \eta$, h为水深, b(x, η)为断面宽度。
根据PSM原理, 在密闭管道上方加一窄缝(见图 1), 有压流也可使用圣维南方程进行计算。当水位高于管顶高程时, 可视其为自由表面流动, 其重力波波速a为
图 1 PSM断面示意
Figure 1. Schematic diagram of cross section in the PSM
$$ a = \sqrt {gA/{T_{\rm{s}}}} $$ (2) 式中: Ts为窄缝宽度。
水击波波速c根据流体特性以及管道材料设定, 通过选择合适的Ts可保证a=c。
$$ {T_{\rm{s}}} = \frac{{gA}}{{{c^2}}} $$ (3) -
采用Godunov格式的有限体积法对控制方程离散, 表达式如下：
$$ \mathit{\boldsymbol{U}}_i^{n + 1} = \mathit{\boldsymbol{U}}_i^n - \frac{{\Delta t}}{{\Delta {x_i}}}({\mathit{\boldsymbol{F}}_{i + 1/2}} - {\mathit{\boldsymbol{F}}_{i - 1/2}}) + \Delta t{\mathit{\boldsymbol{S}}_i} $$ (4) 式中: i为计算单元序号；n为时刻；Δxi为单元i的长度；Δt为计算时间步长, Δt≤Δx/Sw, maxn, Sw, maxn为n时刻整个计算域内最大波速；Fi+1/2和Fi-1/2分别为i+1/2和i-1/2界面处的通量。
-
圣维南方程的Riemann问题可定义为分段变量恒定的初值问题(Initial Value Problem, IVP), 当单元界面处初始条件为一个间断时, 即构成了经典的Riemann问题(Riemann Problem, RP)。表达式如下：
$$ \frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{U}})}}{{\partial x}} = 0\quad \mathit{\boldsymbol{U}}(x,0) = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{U}}_{\rm{L}}}}&{x < 0}\\ {{\mathit{\boldsymbol{U}}_{\rm{R}}}}&{x > 0} \end{array}} \right. $$ (5) 式中: L和R分别代表左和右。
IVP问题可获得精确解并可在x-t坐标系下表示出来(见图 2)。两侧的波一个向左传播, 另一个向右传播, 可能是激波或稀疏波。左波和右波可将坐标系划分为3个区域, 分别为UL、U*和UR, 其中U*为左右不同状态间衍生出的新状态。本文采用Kerger等[3]推导的圣维南方程精确Riemann求解器求解。
图 2 Riemann问题波的结构
Figure 2. Wave configuration in the solution of the Riemann problem
A*为衍生的新状态中过流断面面积, A*的求解是精确Riemann求解器的核心, A*满足如下Riemann问题代数恒等式:
$$ f({A_*}) \equiv {f_{\rm{L}}}({A_*},{A_{\rm{L}}}) + {f_{\rm{R}}}({A_*},{A_{\rm{R}}}) + {u_{\rm{R}}} - {u_{\rm{L}}} = 0 $$ (6) 式中: u为流速；fL和fR为波的界面通量, 表达式如下：
$$ \begin{array}{l} {f_{\rm{L}}} = \left\{ {\begin{array}{*{20}{l}} {\phi ({A_*}) - \phi ({A_{\rm{L}}})}&{{A_*} \le {A_{\rm{L}}}({\rm{稀疏波}})}\\ {\sqrt {\frac{{[g{I_1}({A_*}) - g{I_1}({A_{\rm{L}}})]({A_*} - {A_{\rm{L}}})}}{{{A_*}{A_{\rm{L}}}}}} }&{{A_*} > {A_{\rm{L}}}({\rm{激波}})} \end{array}} \right.\\ {f_{\rm{R}}} = \left\{ {\begin{array}{*{20}{l}} {\phi ({A_*}) - \phi ({A_{\rm{R}}})}&{{A_*} \le {A_{\rm{R}}}({\rm{稀疏波}})}\\ {\sqrt {\frac{{[g{I_1}({A_*}) - g{I_1}({A_{\rm{R}}})]({A_*} - {A_{\rm{R}}})}}{{{A_*}{A_{\rm{R}}}}}} }&{{A_*} > {A_{\rm{R}}}({\rm{激波}})} \end{array}} \right. \end{array} $$ (7) 式中:$\phi(A)=\int_{0}^{A} \sqrt{g /(\alpha b(\alpha))} \mathrm{d} \alpha$, b(α)为断面宽度函数。
Q*为衍生新状态中的流量, Q*的表达式为
$$ {Q_*} = \frac{{{A_*}}}{2}({u_{\rm{L}}} + {u_{\rm{R}}}) + \frac{{{A_*}}}{2}[{f_{\rm{R}}}({A_*},{A_{\rm{R}}}) - {f_{\rm{L}}}({A_*},{A_{\rm{L}}})] $$ (8) -
式(6)是关于A*的非线性方程, 通过迭代可获得其近似解。文献[3]和[12]中均采用Newton-Raphson法迭代求解, 该算法具备二阶收敛, 但缺点是要对fL和fR求偏导, 断面几何形状规则情况下可求出f′L和f′R的表达式, 但复杂几何条件下对fL和fR求偏导几乎不可能。另外, Newton-Raphson法要求函数充分光滑, 而在PSM中, 式(6)在Amax处连续但不可导, 该点为一尖点, 导致函数不光滑(见图 3), Newton-Raphson法往往无法收敛至真实解, 在以往的研究中忽略了这一点, 由此也可导致数值振荡。
图 3 PSM中I1(A)和ϕ(A)示意
Figure 3. Schematic diagram of I1(A) and ϕ(A) in the PSM
本文采用三阶收敛迭代方法结合二分法对式(6)进行迭代求解。裕静静和江平[13]提出的免求导、三阶收敛迭代方法表达式如下：
$$ \left\{ {\begin{array}{*{20}{l}} {{y^n} = A_*^n - \frac{{2{f^2}(A_*^n)}}{{f(A_*^n + f(A_*^n)) - f(A_*^n - f(A_*^n))}}}\\ {A_*^{n + 1} = A_*^n - \frac{{f(A_*^n)}}{{f({y^n}) - f(A_*^n)}}({y^n} - A_*^n)} \end{array}} \right. $$ (9) 当满足如下条件迭代停止
$$ \Delta A = \frac{{|{A^{n + 1}} - {A^n}|}}{{({A^{n + 1}} + {A^n})/2}} < {\varepsilon _1} $$ (10) 式中: ε1为设定的误差容许值。
当迭代次数达到设定值时仍不收敛, 则根据三阶收敛迭代方法确定的上界A*up和下界A*down, 采用二分法继续迭代, 将高阶收敛迭代方法与尖点附近低阶收敛迭代方法相结合, 二分法表达式如下：
$$ A_*^i = (A_*^{{\rm{ up }}} + A_*^{{\rm{ down }}})/2,\left\{ {\begin{array}{*{20}{l}} {A_*^{{\rm{ up }}} = A_*^i}&{f(A_*^{{\rm{ up }}})f(A_*^i) \ge 0}\\ {A_*^{{\rm{ down }}} = A_*^i}&{f(A_*^{{\rm{ down }}})f(A_*^i) \ge 0} \end{array}} \right. $$ (11) 当|A*up-A*down|＜ε2时, 则迭代停止, ε2为设定误差容许值。
-
以上求解步骤可获得A*和Q*, 然而Godunov格式的变量随时间更新需要获得沿x/t=0轴(单元界面处)Riemann问题的解, 具体参见Toro[12]的抽样方法。
采用二阶龙格-库塔离散的时间分裂方法来处理源项[3]：
$$ \mathit{\boldsymbol{U}}_i^{{\rm{adv}}} = \mathit{\boldsymbol{U}}_i^n - \frac{{\Delta t}}{{\Delta x}}[{\mathit{\boldsymbol{F}}_{i + 1/2}} - {\mathit{\boldsymbol{F}}_{i - 1/2}}],\mathit{\boldsymbol{U}}_i^{{\rm{int}}} = \mathit{\boldsymbol{U}}_i^{{\rm{adv}}} + \frac{{\Delta t}}{2}\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{U}}_i^{{\rm{adv}}}),\mathit{\boldsymbol{U}}_i^{n + 1} = \mathit{\boldsymbol{U}}_i^{{\rm{adv}}} + \Delta t\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{U}}_i^{{\rm{int}}}) $$ (12) 有两种类型的边界条件：①固壁边界, 引入反射边界条件, 边界处的水位与邻近单元一致, 而流速与邻近单元相反；②允许波穿过的边界条件, 边界处的水位和流速与邻近单元一致。
-
相关研究表明, 在激波间断与下游状态衔接上, 守恒方程强制以直线衔接, 而零波强方程则以曲线衔接, 这种矛盾是激波附近产生非物理振荡的根本原因[8]。
以充管水流运动为例(见图 4)。在守恒条件下, i单元水头hi和单宽流量qi满足如下关系式[8]：
图 4 管道充水波运动及状态连接示意
Figure 4. Schematic diagram of pipe filling wave movement and state connection
$$ {q_i} = {q_{i - 1}} + {S_{\rm{w}}}({h_i} - {h_{i - 1}}) $$ (13) 式中: qi-1为i-1单元单宽流量；Sw为波速；hi-1为i-1单元水深。
在激波前后零波强条件下(图 4(f)), i单元水头hi和单宽流量qi满足如下关系式[8]：
$$ {q_i} = \frac{{\sqrt {{h_i}} }}{{2{h_{i - 1}}}}2{q_{i - 1}}\sqrt {{h_i}} + ({h_i} - {h_{i - 1}})\sqrt {2g{h_{i - 1}}\frac{{h_i^2 - 2{h_{i - 1}}Z - h_{i - 1}^2}}{{{h_i} - {h_{i - 1}}}}} $$ (14) 式中: Z为超压水头。
在零波强条件下(图 4(f)), qi和hi呈非线性关系, 这与守恒条件下的线性关系存在较大差别。进一步分析, Godunov格式下单元i的n+1时刻变量值为单元内的平均值, 在强间断条件下, 这种变量取平均会产生较大误差。
以T=T2时刻为例(图 4(b)), 图中蓝色为实际物理状态, 而虚线是计算得出的虚假物理状态, 若采用一阶精度变量空间重构方法, 即UiL=UiR=Ui(UiL和UiR分别为i单元左侧和右侧变量重构值), 实际的物理状态下, i+1/2界面处流速为0, 即qi+1/2=0, hi+1/2=hi+1；但经过式(4)计算后, qi+1/2 >0, hi+1/2>hi+1；i-1/2界面处实际物理状态与式(4)计算结果亦有差别。
-
在Godunov格式一阶精度变量空间重构方法下, 非物理振荡因激波间断在单元格内被人为拉平所致, 不能反映激波间断下水流运动的真实物理状态。高精度格式亦不能很好地抑制非物理振荡, 加大数值黏性抑制振荡方法人为因素大。对变量进行精确的空间重构是抑制非物理振荡的关键, 而精确Riemann解可准确计算变量在单元内的空间分布状态, 因此, 本文提出基于精确Riemann解的变量空间重构方法。
明满流过渡模拟中的非物理振荡主要产自于激波前缘所在的计算单元, 将该单元的一阶精度变量重构方法改为基于精确Riemann解的变量空间重构方法, 其他单元不变。设激波前缘所在单元i, 当前时刻为n, 经过计算时间步长Δt后的时刻为n+1, 当Δt≤Δx/Sw, maxn, 在Δt内根据抽样方法, 存在如下表达式[14]：
$$ \mathit{\boldsymbol{\tilde U}}({x_{i - 1/2}},t) = {\mathit{\boldsymbol{U}}_{i - 1/2}}(0) = {\rm{ constant }}\quad \mathit{\boldsymbol{\tilde U}}({x_{i + 1/2}},t) = {\mathit{\boldsymbol{U}}_{i + 1/2}}(0) = {\rm{constant}} $$ (15) 式中: $\mathit{\boldsymbol{\tilde U}}\left( {{x_{i + 1/2}}, t} \right)$为i+1/2处变量随时间变化函数；Ui+1/2(0)为Riemann问题RP(Uin, Ui+1n)沿x/t =0轴的解。
根据式(15)可得如下推论：
$$ \mathit{\boldsymbol{U}}_{i - 1/2}^n = \mathit{\boldsymbol{U}}_{i - 1/2}^{n + 1} = {\mathit{\boldsymbol{U}}_{i - 1/2}}(0)\quad \mathit{\boldsymbol{U}}_{i + 1/2}^n = \mathit{\boldsymbol{U}}_{i + 1/2}^{n + 1} = {\mathit{\boldsymbol{U}}_{i + 1/2}}(0) $$ (16) i单元n+1时刻变量空间重构表达式如下:
$$ \mathit{\boldsymbol{U}}_i^{\rm{L}} = \mathit{\boldsymbol{U}}_{i - 1/2}^{n + 1},\mathit{\boldsymbol{U}}_i^{\rm{R}} = \mathit{\boldsymbol{U}}_{i + 1/2}^{n + 1}\quad A_i^{n + 1} \le \max (A_i^{\rm{L}},A_i^{\rm{R}}) $$ (17) 式中: Ain+1为i单元n+1时刻平均过流断面面积；AiL和AiR分别为i单元左侧和右侧过流断面面积重构值。根据精确Riemann解确定的i单元变量重构结果, 仍然采用式(4)对i单元的变量在时间维上更新, 也就意味着变量空间重构此时不依赖Uin+1的计算结果, 避免了变量单元平均对过渡过程计算精度的损害。
当Ain+1>max(AiL, AiR)时, 表明激波前缘在Δt内已经通过i+1/2界面而进入i+1单元, 这时需要“破时段”, 通过试算确定Δt′, 即经过Δt′, $\left|A_{i}^{n+\Delta t^{\prime} / \Delta t}-\max \left(A_{i}^{\mathrm{L}}, A_{i}^{\mathrm{R}}\right)\right|<\varepsilon_{3}$, ε3为容许误差, 间断前缘位置与i+1/2界面重合；在此基础上计算经Δt″(Δt″=Δt-Δt′)后Uin+1和Ui+1n+1的结果。
-
为了便于将数值计算结果与理论值进行对比, 给出一种理想化试验装置[6]。400 m长光滑、水平放置矩形试验装置, 断面尺寸为1 m×1 m, 试验装置内初始水位(水深)为0.6 m, 初始流速为0 m/s。试验装置上、下游均与水库连接, 下游水库水位恒定, 为0.6 m, 上游水库水位从初始0.6 m突然增高至4 m, 试验装置示意见图 5, 设水击波波速为100 m/s(与文献[6]保持一致)。
图 5 充水试验装置示意
Figure 5. Schematic diagram of pipe filling experimental apparatus
上游水库和管道入口处满足能量守恒方程, 间断处波速按照右侧激波波速计算, 表达式如下：
$$ {{H_{{\rm{res}}}} = {H_{{\rm{pipe}}}} + \frac{{u_{{\rm{pipe}}}^2}}{{2g}}} $$ (18) $$ {{S_{\rm{w}}} = {u_{\rm{R}}} + \sqrt {\frac{{[g{I_1}({A_*}) - g{I_1}({A_{\rm{R}}})]{A_*}}}{{({A_*} - {A_{\rm{R}}}){A_{\rm{R}}}}}} } $$ (19) 式中: Hres为上游水库水位；Hpipe为管道入口处水头；upipe为管道入口处平均流速；L、R表示间断处左侧和右侧。
在有压流下, 管道水头的理论值为3.167 m, 流速为4.044 m/s, 波速的理论值为10.077 m/s。图 6给出了一阶精度空间重构方法下x=20 m处计算水头随时间变化过程。理论上当t=1.98 s时, 间断前缘到达x=20 m处, 此后该处产生比较明显的非物理振荡：单元长度为4 m时, 振荡幅度较大, t>18.5 s后振荡消失, 计算结果收敛于理论值；单元长度为1 m时, 振荡幅度较小, t>3.2 s后振荡消失, 计算结果收敛于理论值。
图 6 x =20 m处管道水位计算值与理论值对比(一阶精度空间重构法)
Figure 6. Comparison of calculated and theoretical water level of pipeline at x =20 m(first order precision)
图 7给出了30 s时基于一阶精度变量空间重构方法的计算结果, 并与理论值以及文献[6]的计算结果进行对比。从图中可见, 本文波速、水头、流速的模拟结果均较准确, 收敛速度快, 非物理振荡只在激波间断附近产生, 并且随着计算单元尺寸的减小, 振幅缩小, 计算结果逐渐趋于理论值；文献[6]与图 7(a)的计算条件完全一致, 其计算结果收敛速度慢, 靠近间断前缘处数值振荡较强烈, 水头最大值达到了4.95 m, 振荡幅度远高于本文计算结果；二者计算结果的差异可能是因为文献[6]在采用Newton-Raphson法对Riemann问题代数恒等式迭代求解过程中, 忽略了尖点(A =Amax=1 m2)存在而导致迭代未收敛至真实解所造成。
图 7 30 s时管道水位计算值与理论值对比(一阶精度空间重构法)
Figure 7. Comparison of calculated and theoretical water level of pipeline in 30 s(first order precision)
图 8给出了30 s时基于精确Riemann解的变量空间重构方法的计算结果。从图中可见, 计算结果与理论值较为一致, 从机理上消除了计算过程中的非物理振荡。在Godunov格式中, 对变量进行空间精确重构, 可保证单元界面处通量计算的准确性, 进而可保证计算结果的精度。本文提出的基于精确Riemann解的变量空间重构方法实现了变量空间精确重构, 采用三阶收敛算法结合二分法对Riemann问题代数恒等式进行迭代求解, 避免了因尖点存在导致无法收敛至真实解, 这些都是抑制非物理振荡的关键。
图 8 30 s时管道水位计算值与理论值对比(基于精确Riemann解空间重构法)
Figure 8. Comparison of calculated and theoretical water level of pipeline in 30 s(based on exact Riemann solution)
-
Wiggert试验装置为一水平长30 m、宽0.51 m的水槽, 在水槽中间安装10 m长的顶盖, 使其形成一个0.51 m宽、0.148 m高的密闭矩形管道[3, 15], 见图 9。水槽床面为光滑混凝土, 参照文献[3]和[15], 糙率n取0.01, 水击波波速取20 m/s, 初始条件为水位0.128 m的静水。左侧波使水槽带顶盖部位形成有压管流, 在水槽装置中布置了4个压力测点(a、b、c、d)。Wiggert在顶盖上下游各布置了1个水位观测点, 本研究将上下游水位随时间变化过程作为数值计算的边界条件(图 10(a), 图 10(b))。
图 10 Wiggert试验边界条件及计算结果
Figure 10. Wiggert test boundary condition and results
从图 10中可见, 上游端2 s左右开始明满流过渡, 过渡过程在沿管道向下游发展过程中, 间断波前缘越来越陡, 于7 s左右到达下游断面。4个测点计算结果与实测值均比较吻合, 间断波前缘位置与实测值比较接近, 模拟结果明显优于相同试验下TPA模型模拟结果[15]。TPA模型采用Roe近似Riemann求解器, 将非线性问题线性化求解, 另外对强间断所在单元未作抑制数值振荡的技术处理。
-
(1) 采用Preissmann窄缝法对明满流过渡过程进行模拟, 基于Godunov格式, 采用精确Riemann求解器求解, 该求解器是一种非线性、全波族的求解器, 对于明满流过渡过程模拟具有较强的适用性。
(2) Preissmann窄缝法中, Riemann问题代数恒等式在明满流交界处连续但不可导, 存在尖点, 采用Newton-Raphson及类似加速收敛方法均不能保证收敛至真实解。本研究采用免求导、三阶收敛算法结合二分法进行迭代求解, 可稳定收敛至真实解, 特别适用于不光滑函数条件下非线性方程根的迭代求解。但该方法收敛速度较慢, 如何提高迭代收敛速度将是下一步研究的重点。
(3) 非物理振荡是限制Preissmann窄缝法应用的关键因素, 因过渡过程存在较强的激波间断, 现有的变量空间重构方法很难准确表达真实的物理状态, 导致数值解出现非物理振荡。本文提出了基于精确Riemann解的变量空间重构方法, 准确表达强间断水流在单元内的空间分布状态, 保证了单元界面处通量计算的精度, 从机理上抑制了非物理振荡。
(4) 实例研究表明, 计算值与解析解或实测值吻合良好, 研究成果为明满流过渡过程的高精度数值模拟提供了新的方法。
Simulation of transient mixed flows based on the exact Riemann solver
-
摘要: Preissmann窄缝法模拟明满流过渡过程方法简单，但存在明显的非物理振荡，抑制非物理振荡是该方法应用的关键。基于Godunov格式和精确Riemann求解器对明满流过渡过程进行模拟，针对Riemann问题代数恒等式在明满流交界处不光滑问题，提出了三阶收敛方法与二分法结合的迭代求解方法，保证迭代收敛至真实解；针对由于变量空间重构方法不能准确表达变量在空间中真实物理状态而导致的非物理振荡，提出了基于精确Riemann解的变量空间重构方法，准确表达激波间断在单元内的空间分布状态，从机理上抑制了非物理振荡。实例研究表明，数值计算结果与解析解或实测值吻合良好，研究成果为明满流过渡过程的高精度数值模拟提供了新的方法。
-
关键词:
- 明满流过渡 /
- 非物理振荡 /
- Preissmann窄缝法 /
- Godunov格式 /
- 精确Riemann求解器
Abstract: The Preissmann slot method is simple to simulate transient mixed flows, but it has obvious non-physical oscillation. The key to the application of this method is to suppress the non-physical oscillations. The transient mixed flows are simulated based on the Godunov scheme and the exact Riemann solver. In order to solve the problem that the algebraic identity of the Riemann problem is not smooth at the boundary of open and pressure flow, an iterative solution method is proposed. This method combines the third-order convergence method and the dichotomy method to ensure the iterative convergence to the real solution. Considering the fact that the non-physical oscillation caused by the variable space reconstruction method cannot accurately express the real physical state of the variable in space, a variable space reconstruction method based on the exact Riemann solution is proposed. This method can accurately express the spatial distribution state of the shock discontinuity in the cell, and the non-physical oscillation is suppressed from the mechanism. The case study shows that the numerical results are in good agreement with the analytical solutions or the measured values. The research results provide a new method for the high-precision numerical simulation of the transient mixed flows. -
图 4 管道充水波运动及状态连接示意
Figure 4. Schematic diagram of pipe filling wave movement and state connection
图 6 x =20 m处管道水位计算值与理论值对比(一阶精度空间重构法)
Figure 6. Comparison of calculated and theoretical water level of pipeline at x =20 m(first order precision)
图 7 30 s时管道水位计算值与理论值对比(一阶精度空间重构法)
Figure 7. Comparison of calculated and theoretical water level of pipeline in 30 s(first order precision)
图 8 30 s时管道水位计算值与理论值对比(基于精确Riemann解空间重构法)
Figure 8. Comparison of calculated and theoretical water level of pipeline in 30 s(based on exact Riemann solution)
-
[1] PACHALY R L, VASCONCELOS J G, ALLASIA D G, et al. Comparing SWMM 5.1 calculation alternatives to represent unsteady stormwater sewer flows[J]. Journal of Hydraulic Engineering, 2020, 146(7):04020046. doi: 10.1061/(ASCE)HY.1943-7900.0001762 [2] VASCONCELOS J G, WRIGHT S J, ROE P L. Two-component pressure approach for the simulation of flow regime transition in sewers[J]. Journal of Hydraulic Engineering, 2006, 132(6):553-562. doi: 10.1061/(ASCE)0733-9429(2006)132:6(553) [3] KERGER F, ARCHAMBEAU P, ERPICUM S, et al. An exact Riemann solver and a Godunov scheme for simulating highly transient mixed flows[J]. Journal of Computational and Applied Mathematics, 2011, 235(8):2030-2040. doi: 10.1016/j.cam.2010.09.026 [4] 王志刚.满足3个守恒律的Godunov型格式[J].上海交通大学学报, 2012, 46(10):1697-1706. WANG Z G. Godunov type scheme satisfying three conservation laws[J]. Journal of Shanghai Jiao Tong University, 2012, 46(10):1697-1706. (in Chinese) [5] LEÓN A S, GHIDAOUI M S, SCHMIDT A R, et al. Application of Godunov-type schemes to transient mixed flows[J]. Journal of Hydraulic Research, 2009, 47(2):147-156. doi: 10.3826/jhr.2009.3157 [6] MALEKPOUR A, KARNEY B W. Spurious numerical oscillations in the Preissmann slot method:origin and suppression[J]. Journal of Hydraulic Engineering, 2016, 142(3):04015060. doi: 10.1061/(ASCE)HY.1943-7900.0001106 [7] MAO Z H, GUAN G H, YANG Z H. Suppress numerical oscillations in transient mixed flow simulations with a modified HLL solver[J]. Water, 2020, 12(5):1245. doi: 10.3390/w12051245 [8] VASCONCELOS J G, WRIGHT S J, ROE P L. Numerical oscillations in pipe-filling bore predictions by shock-capturing models[J]. Journal of Hydraulic Engineering, 2009, 135(4):296-305. doi: 10.1061/(ASCE)0733-9429(2009)135:4(296) [9] AN H, LEE S, NOH S, et al. Hybrid numerical scheme of Preissmann slot model for transient mixed flows[J]. Water, 2018, 10(7):899. doi: 10.3390/w10070899 [10] TORO E F, GARCIA-NAVARRO P. Godunov-type methods for free-surface shallow flows:a review[J]. Journal of Hydraulic Research, 2007, 45(6):736-751. doi: 10.1080/00221686.2007.9521812 [11] CHAUDHRY M H. Open-channel flow[M]. Boston:Springer, 2008:333-339. [12] TORO E F. Shock-capturing methods for free-surface shallow flows[M]. Chichester:John Wiley & Sons, 2001:93-108. [13] 裕静静, 江平.不用求导含参数的三阶收敛迭代方法[J].大学数学, 2015, 31(3):34-38. doi: 10.3969/j.issn.1672-1454.2015.03.007 YU J J, JIANG P. Third-order convergent iterative method with parameter without derivation[J]. College Mathematics, 2015, 31(3):34-38. (in Chinese) doi: 10.3969/j.issn.1672-1454.2015.03.007 [14] TORO E F. Riemann solvers and numerical methods for fluid dynamics[M]. Heidelberg:Springer Berlin Heidelberg, 2009:213-218. [15] VASCONCELOS J G, WRIGHT S J. Comparison between the two-component pressure approach and current transient flow solvers[J]. Journal of Hydraulic Research, 2007, 45(2):178-187. doi: 10.1080/00221686.2007.9521758 -
点击查看大图
图(10)
计量
- 文章访问数: 13
- HTML全文浏览量: 2
- PDF下载量: 6
- 被引次数: 0
