一种求解浅水问题模拟间断水流数值的方法
【专利摘要】本发明公开了一种求解浅水问题模拟间断水流数值的方法，涉及计算流体力学领域，具体公开了一种求解二维双曲型方程组黎曼问题的数值方法，将其应用于求解浅水波方程，进行间断水流的数值模拟。本发明在二维无结构三角形网格中，求解双曲型浅水方程的黎曼问题，提高了数值精度和分辨率，探讨了间断水流的水动力特征。本发明基于无结构三角形网格建立了能够适应复杂地形和复杂计算区域中带有大梯度或间断水流现象的高分辨率数值模型，该方法可以较大程度的减小误差，提高模拟水动力特征的分辨率。
【专利说明】一种求解浅水问题模拟间断水流数值的方法
【技术领域】
[0001] 本发明公开了一种求解浅水问题模拟间断水流数值的方法，具体涉及计算流体力 学【技术领域】。
【背景技术】
[0002] 数值模拟是了解水动力特征的主要途径之一。在近海、河口、湖泊、大型水库等大 面积开敞水域，可以采用二维浅水波方程来描述水动力参数的特征。由于浅水波方程具有 非线性特性，一般情况下无法获得解析解，主要通过数值模拟求解。
[0003] 数值方法是数值模拟的核心。河口、海岸及近海水域，受外海潮波、岸线和地形等 诸多因素的影响，水动力环境十分复杂，除了这些物理因素外，还有浅水波方程本身的数学 理论特性，如果数值方法选择不当会导致虚假震荡，并且不是所有适合求解浅水水流的数 值方法都可以用来模拟间断水流。一般来说，即使给定的初值条件非常简单，浅水波方程的 解也可能发展成间断波或激波，事实上，这是无耗散性的非线性双曲守恒率的一般特性。
[0004] 溃坝波、涌潮等强间断流动的数值模拟始终是水动力学研究中重要的、也是困难 的课题。很多传统的求解浅水波方程的数值方法都是基于微分形式的控制方程，要求物理 变量处处连续、可微，因此在模拟这类水位、流速存在间断的流动时便会出现一些问题。国 内外一些通用水动力学数值模拟软件，应用于一般河道都很成功，但是还不能很好的模拟 间断流动，特别是强间断流动。间断流动数值模拟的主要困难是建立"和谐"的计算格式， 即浅水方程通量项要与源项"和谐"；间断前后水位、流速梯度很大，要求计算格式既能模拟 大梯度流动，又要格式稳定。因而，间断解的理论研究和数值模拟就具有特别重要的学术意 义。
[0005] 理论分析表明，利用高阶、高精度数值格式不仅可以降低对网格的要求，还能够准 确捕捉流场中复杂的流动现象，因此发展高精度数值方法是提高数值模拟计算精度的重要 途径之一。另外，对于空间和时间存在多种尺度变化的物理现象的数值模拟，要求多种尺度 信息在数值计算中被捕捉到，这也需要高精度高分辨率格式。因此浅水水流中间断流动的 高精度数值模拟更具有重要的学术价值。
[0006] 正确的模拟间断水流首先需要正确的模拟间断关系。在不考虑粘性效应时，间断 波是一间断面。当前较有效的方法是，把间断面包含在计算区域内用统一的方法计算，即所 谓的激波捕捉法。其基本出发点是使用与守恒律微分方程组相容的守恒性差分格式，所得 的差分解在间断两侧自动满足间断条件，因而不论解中是否存在间断，可以不加区别地统 一进行计算，但要求计算网格较密，且要求计算格式具有模拟大梯度的能力。常见的捕捉法 中的一阶精度格式，会使激波抹的过平，二阶或高于二阶精度的经典格式数值解在激波附 近易产生非物理震荡。近年来为了改善激波附近数值解的分辨率和虚假震荡取得了很大进 展，即产生了一些求解浅水间断波的高分辨率数值格式。
[0007] 所谓"高分辨率"指数值解图形的锐利和逼真性，表明数值方法的数值耗散微弱且 适中。以忽略源项的守恒型浅水波方程U t+F(U)x = 0为例，其中U为守恒性物理量，F(U) 为通量项。设真解的间断点xs在区间+.A',, j内，其中下标i为网格节点编号，当x < xs 时U =队，当X > xs时U = UK，带下标L, R的符号Uu UK分别表示U在xs点左边和右边的 值。
[0008] 守恒格式的差分解为： 卜（./<0
[0009] U- - -,UU (7 - /); \UR (./> 0
[0010] 其中，下标j为网格节点编号，仏为\点上的近似U值。
[0011] 为了满足守恒性，需要下式成立： _2] Uy (X +, ^^) = UL (.τ, ^x,H) + UR (x+ ； ^)-
[0013] 此式和差分格式唯一地确定了XJPUM。因此，数值间断必须至少涉及三个格点（即 两个区间）。目前一些优良格式的间断分辨率已达到二格至三格。要实现高分辨率格式的 设计，就是在保证格式的相容性、稳定性的前提下，在设计和构造格式时能够尽可能减弱其 数值耗散性，保持数值耗散性优势。
[0014] 网格生成是数值模拟的基础，它实质是物理求解域与计算求解域的转换。在求解 具有复杂几何形状的流场时，适当的网格生成是一个十分重要的问题，网格质量的好坏直 接影响到计算结果的收敛情况及精度。无结构网格具有复杂区域适应性好、局部加密灵活 和便于自适应的优点，能很好地模拟自然边界及复杂的水下地形，提高边界模拟精度。
【发明内容】
[0015] 本发明所要解决的技术问题是：针对现有技术的缺陷，提供一种求解浅水问题模 拟间断水流数值的方法，，具体是一种求解二维浅水波方程黎曼解的数值方法，可以应用于 模拟溃坝波、涌潮等不可压间断水流，进而了解间断水流的水动力特征。
[0016] 本发明为解决上述技术问题采用以下技术方案：
[0017] 一种求解浅水问题模拟间断水流数值的方法，建立模拟间断水流问题的高分辨率 数值模型，利用所述模型模拟地形和计算域中带有大梯度解或间断的水流现象，具体步骤 如下：
[0018] 步骤一、建立计算网格，将物理计算区域进行无结构三角形网格分割，得到复数个 三角形网格单元，将单一的网格单元作为控制元，其中第i个三角形控制元记为;
[0019] 步骤二、通过守恒型物理量反映每个计算网格的水动力特征，在笛卡尔直角坐标 系下，将平面二维浅水运动守恒形式的方程表示为：
[0020] Ut+WU)x+G^y = S (1);
[0021] 其中，U为守恒量向量；F，G为通量向量；S为源项向量，t为时间，X和y为笛卡尔 坐标系坐标，U t，F(U)X，G(U)y分别表示为时间t和空间X，y方向的偏导数，具体表示为：
[0022] U = [h, uh, vh]T ；
[0023] F = [uh,u2h+gh2/2,uvh]T ；
[0024] G = [vh, uvh, v2h+gh2/2]T ；
[0025] S = [0, -ghbx, -ghby]T ；
[0026] 其中，T表示向量的转置；u、v分别是x、y方向沿水深平均的流速分量，并且为 关于X，y，t的函数；h是水深；g是重力加速度；b x、by分别为河床比降在X、y方向上的分 量，b (X，y)为河底高程函数；
[0027] 求解上述方程的黎曼解，对方程组（1)，令向量函数E= (F，G)，V= 为梯 ^ cx cn- J 度算子，将式（1)表示为：
[0028] Ut+ ▽ · E = S (2)；
[0029] 采用步骤一中划分的三角形网格单元对计算区域进行离散，在上对方程组（2) 进行积分，并按逆时针利用Green公式将面积分化为线积分，得：
【权利要求】
1. 一种求解浅水问题模拟间断水流数值的方法，其特征在于，建立模拟间断水流问题 的高分辨率数值模型，利用所述模型模拟地形和计算域中带有大梯度解或间断的水流现 象，具体步骤如下： 步骤一、建立计算网格，将物理计算区域进行无结构三角形网格分割，得到复数个三角 形网格单元，将单一的网格单元作为控制元，其中第i个三角形控制元记为; 步骤二、通过守恒型物理量反映每个计算网格的水动力特征，在笛卡尔直角坐标系下， 将平面二维浅水运动守恒形式的方程表示为： Ut+F(U)x+G(U)y = S (1); 其中，U为守恒量向量；F，G为通量向量；S为源项向量，t为时间，X和y为笛卡尔坐标 系坐标，Ut，F(U)X，G(U)y分别表示为时间t和空间X，y方向的偏导数，具体表示为 ： U = [h, uh, vh]T ； F = [uh, u2h+gh2/2, uvh]T ； G = [vh, uvh, v2h+gh2/2]T ； S = [0, -ghbx, -ghby]T ； 其中，T表示向量的转置；u、v分别是x、y方向沿水深平均的流速分量，并且为关于 X，y，t的函数；h是水深；g是重力加速度；bx、by分别为河床比降在x、y方向上的分量； 求解上述方程的黎曼解，对公式（1)，令向量函数E = (F，G)，▽ 为梯度算 ^ ox oy ) 子，将式⑴表示为： Ut+ ▽ · E = S (2)； 采用步骤一中划分的三角形网格单元对计算区域进行离散，在上对方程组（2)进 行积分，并按逆时针利用Green公式将面积分化为线积分，得： j 3 UdQ = -Y ? E-nJI+? Sdil (3)； di tfiL·-·· Jn.+ 其中，为时间方向的常微分算子，cm是面积分微元，di是线积分微元，Lm表示第i di 个三角形单元Ω i的第m条边，m = 1，2, 3, nm = (nmx, nmy) = (cos θ，sin θ )为Lm的单位外 法向量，nmx，nmySnm在x、y方向的分量，Θ为外法线方向与水平方向的夹角； 定义单元平均值士由1严，&赫，亂其中，ΙΩ」为网格单元面积； 利用两点Gauss求积公式求解式（3)中的线积分： I ￡ ·njl?|LmImqEψ nm ⑷； m q=i 其中，I Lm I为三角形边界线Lm的长度，G,为三角形单元边Lm上的Gauss点，(f = u) _ 1 为相应的权系数，q为咼斯点编号，?I = 0? ; 构造一个数值通量格式近似替代E(U(Gq，t))，将三角形网格单元边界上内部和外部 的变量值分别记为队、UK，将通过Gauss点G,的通量记为，其中， UjGq，t)，UjGpt)分别为U在Gauss点三角形网格单元内部和外部的值，由此得出控制方 程（3)的有限体积半离散化格式：
步骤三、求解三角形网格单元边界上的HLL数值通量# ; 步骤四、采用三阶Runge-Kutta时间离散方法对时间方向上的离散进行处理，得出三 阶精度的时间离散格式； 步骤五、采用三阶WENO格式重构方法求解UJG,，t)，UJG,，t); 步骤六、对每个网格单元，利用当前时间层上每个单元的守恒量向量积分平均值，求解 下一时刻的守恒量向量积分平均值： 经过步骤五对流体方程时间方向上的的离散后，再对流体方程进行空间离散，表达式 为： ^ = L(U)； di 其中，L(U)是空间离散化算子； 设定第η个时间层上的U值为Un，则经过一个时间步长Λ t后，下个时间层的近似值 Un+1表示为： U(1) = Un+AtL(Un)；
其中，U(1)、U(2)为中间过渡向量。
2.如权利要求1所述的一种求解浅水问题模拟间断水流数值的方法，其特征在于，所 述步骤三中，方程的HLL通量表达式具体如下：
其中，sK分别是波速度的上下确界At = %-a!^，sK = uK+aKqK ; 其中，左右两边的波速表示为：? = =^Γ; 设定qK、f为中间过渡量，Κ = L，R，具体为：
3.如权利要求1所述的一种求解浅水问题模拟间断水流数值的方法，其特征在于，所 述步骤四包括：针对三角形网格单元构造三角形单元的多项式，具体步骤包括： (501) 线性多项式的构造： 基于三角形单元以及相邻的九个三角形单元，设定九个模板，具体为： Bi = {Ω0, Ωρ QJ ；B2 = {Ω0, Qk, QJ ；B3 = {Ω0, Ω?； Ω j} ；B4 = {Ω0, Ω?； Ω?3}； Β5 = {Ω0, Ω?； Qib} ；Β6 = {Ω0, Qj, QJa} ；Β7 = {Ω0, Ω j, Ω Jb} ；Β8 = {Ω0, Ωk, ΩkJ ；Β9 = { Ω 0, Ω k, Ω kb}； 其中，下标i, j, k, ia, ib, ja, jb, ka, kb为不同相邻单元的编号； 定义T = U B2 U…U B9为混合模板集，其中U为并集符号，在每个小模板上构造线 性多项式P(x，y)，使得P(x，y)在小模板所覆盖的三个网格单元上的平均值和已知单元Ω。 的平均值相等； (502) 构造线性权： 构造如下二次多项式：Q(x，y) = h+t^x+t^y+t^xy+by+bjjy2 ; 其中,131，132，133，134，13 5，136为多项式的待定系数，在{〇(|，1，」， 1;，^，如0>0# Qka, Qkb}这十个单元上满足：, s = 0, i, j, k, ia, ib, ja, jb, ka, kb ;其中 ￡/s 表示单元Ω s上U的平均值， 对三角形单元Ω ^的每个Gauss积分点（xG, yG)，设定线性权为γ i (1 = 1，. . .，9), γ 1 是由三角形网格的几何性质唯一决定的常数，权系数非负； 将线性多项式Ρ (X，y)用线性权Y i合并： 伞，) = Σ#(χ，.ν) (7) i=I y i 满足下述等式：R(xe, ye) = Q(xe, ye) (8) (503) 构造非线性权和光滑因子： 在线性权的基础上构造非线性权，基于一致性，权系数和为1，即Σα二】； /~ι 将九个线性多项式分为四组： /=1 /=1 对于i边上的第一个Gauss点： 第一组含有方程 P2 (〇, k，i)，P4 (0, i，ia)，P5 (0, i，ib)， ^ "(/2^2 +^4^4 +/5^5)^(/2 +^4 + 7s) ? f\ =?Ζ+?4+?5' 第二组含有方程 P3(〇, i，j)，P6(〇, j，ja)，P7(0, j，jb)， h = (?3ρ3 + r,A+￥1^)1 (h + r(,+7t)? ?2= n + +r7： 第三组含有方程 Pi (〇, j，k)，P8 (0, k，ka)，P9 (0, k，kb)， k = (--+?Α+υ?)/(λ +λ+λ)> η=7ι+η+ηι 第四组含有方程 Pjo, j，k)，Ρ2(0, k，i)，Ρ3(0, i，j)， hΗ--+?Ρ2+τ?)^{? +Ti +h)* Λ = +r,： 得出线性多项式：坏^卜工方巧卜^^为多项式的非负解； /-I Π 令：) = ^w。+ a2u, + 〇義+α4-- + a5?ift，对一个n次多项式p (X，y)，定义光滑
测量：ψ= Σ ]>Γ'(dXu))2 Ι<--/|<? 其中，a是一个复指数，D是偏导数，则非线性权为
其中，-?是第?个多项式^的光滑度，ε是小正数。
4. 如权利要求3所述的一种求解浅水问题模拟间断水流数值的方法，其特征在于：所 述ε的取值范围为（1(Γ6, ΚΓ2)。
5. 如权利要求4所述的一种求解浅水问题模拟间断水流数值的方法，其特征在于：所 述ε的取值为ε = ΚΓ6。
6. 如权利要求1所述的一种求解浅水问题模拟间断水流数值的方法，其特征在于：所 述步骤六中，时间步长At的具体计算方式为：
其中，CFL表示CFL条件数。
【文档编号】G06F19/00GK104091065SQ201410315720
【公开日】2014年10月8日 申请日期:2014年7月3日 优先权日:2014年7月3日
【发明者】卢长娜 申请人:南京信息工程大学