一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法
【專利摘要】本發(fā)明公開(kāi)了一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，涉及計(jì)算流體力學(xué)領(lǐng)域，具體公開(kāi)了一種求解二維雙曲型方程組黎曼問(wèn)題的數(shù)值方法，將其應(yīng)用于求解淺水波方程，進(jìn)行間斷水流的數(shù)值模擬。本發(fā)明在二維無(wú)結(jié)構(gòu)三角形網(wǎng)格中，求解雙曲型淺水方程的黎曼問(wèn)題，提高了數(shù)值精度和分辨率，探討了間斷水流的水動(dòng)力特征。本發(fā)明基于無(wú)結(jié)構(gòu)三角形網(wǎng)格建立了能夠適應(yīng)復(fù)雜地形和復(fù)雜計(jì)算區(qū)域中帶有大梯度或間斷水流現(xiàn)象的高分辨率數(shù)值模型，該方法可以較大程度的減小誤差，提高模擬水動(dòng)力特征的分辨率。
【專利說(shuō)明】一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法

【技術(shù)領(lǐng)域】
[0001] 本發(fā)明公開(kāi)了一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，具體涉及計(jì)算流體力 學(xué)【技術(shù)領(lǐng)域】。

【背景技術(shù)】
[0002] 數(shù)值模擬是了解水動(dòng)力特征的主要途徑之一。在近海、河口、湖泊、大型水庫(kù)等大 面積開(kāi)敞水域，可以采用二維淺水波方程來(lái)描述水動(dòng)力參數(shù)的特征。由于淺水波方程具有 非線性特性，一般情況下無(wú)法獲得解析解，主要通過(guò)數(shù)值模擬求解。
[0003] 數(shù)值方法是數(shù)值模擬的核心。河口、海岸及近海水域，受外海潮波、岸線和地形等 諸多因素的影響，水動(dòng)力環(huán)境十分復(fù)雜，除了這些物理因素外，還有淺水波方程本身的數(shù)學(xué) 理論特性，如果數(shù)值方法選擇不當(dāng)會(huì)導(dǎo)致虛假震蕩，并且不是所有適合求解淺水水流的數(shù) 值方法都可以用來(lái)模擬間斷水流。一般來(lái)說(shuō)，即使給定的初值條件非常簡(jiǎn)單，淺水波方程的 解也可能發(fā)展成間斷波或激波，事實(shí)上，這是無(wú)耗散性的非線性雙曲守恒率的一般特性。
[0004] 潰壩波、涌潮等強(qiáng)間斷流動(dòng)的數(shù)值模擬始終是水動(dòng)力學(xué)研究中重要的、也是困難 的課題。很多傳統(tǒng)的求解淺水波方程的數(shù)值方法都是基于微分形式的控制方程，要求物理 變量處處連續(xù)、可微，因此在模擬這類水位、流速存在間斷的流動(dòng)時(shí)便會(huì)出現(xiàn)一些問(wèn)題。國(guó) 內(nèi)外一些通用水動(dòng)力學(xué)數(shù)值模擬軟件，應(yīng)用于一般河道都很成功，但是還不能很好的模擬 間斷流動(dòng)，特別是強(qiáng)間斷流動(dòng)。間斷流動(dòng)數(shù)值模擬的主要困難是建立"和諧"的計(jì)算格式， 即淺水方程通量項(xiàng)要與源項(xiàng)"和諧"；間斷前后水位、流速梯度很大，要求計(jì)算格式既能模擬 大梯度流動(dòng)，又要格式穩(wěn)定。因而，間斷解的理論研究和數(shù)值模擬就具有特別重要的學(xué)術(shù)意 義。
[0005] 理論分析表明，利用高階、高精度數(shù)值格式不僅可以降低對(duì)網(wǎng)格的要求，還能夠準(zhǔn) 確捕捉流場(chǎng)中復(fù)雜的流動(dòng)現(xiàn)象，因此發(fā)展高精度數(shù)值方法是提高數(shù)值模擬計(jì)算精度的重要 途徑之一。另外，對(duì)于空間和時(shí)間存在多種尺度變化的物理現(xiàn)象的數(shù)值模擬，要求多種尺度 信息在數(shù)值計(jì)算中被捕捉到，這也需要高精度高分辨率格式。因此淺水水流中間斷流動(dòng)的 高精度數(shù)值模擬更具有重要的學(xué)術(shù)價(jià)值。
[0006] 正確的模擬間斷水流首先需要正確的模擬間斷關(guān)系。在不考慮粘性效應(yīng)時(shí)，間斷 波是一間斷面。當(dāng)前較有效的方法是，把間斷面包含在計(jì)算區(qū)域內(nèi)用統(tǒng)一的方法計(jì)算，即所 謂的激波捕捉法。其基本出發(fā)點(diǎn)是使用與守恒律微分方程組相容的守恒性差分格式，所得 的差分解在間斷兩側(cè)自動(dòng)滿足間斷條件，因而不論解中是否存在間斷，可以不加區(qū)別地統(tǒng) 一進(jìn)行計(jì)算，但要求計(jì)算網(wǎng)格較密，且要求計(jì)算格式具有模擬大梯度的能力。常見(jiàn)的捕捉法 中的一階精度格式，會(huì)使激波抹的過(guò)平，二階或高于二階精度的經(jīng)典格式數(shù)值解在激波附 近易產(chǎn)生非物理震蕩。近年來(lái)為了改善激波附近數(shù)值解的分辨率和虛假震蕩取得了很大進(jìn) 展，即產(chǎn)生了一些求解淺水間斷波的高分辨率數(shù)值格式。
[0007] 所謂"高分辨率"指數(shù)值解圖形的銳利和逼真性，表明數(shù)值方法的數(shù)值耗散微弱且 適中。以忽略源項(xiàng)的守恒型淺水波方程U t+F(U)x = 0為例，其中U為守恒性物理量，F(xiàn)(U) 為通量項(xiàng)。設(shè)真解的間斷點(diǎn)xs在區(qū)間+.A',, j內(nèi)，其中下標(biāo)i為網(wǎng)格節(jié)點(diǎn)編號(hào)，當(dāng)x < xs 時(shí)U =隊(duì)，當(dāng)X > xs時(shí)U = UK，帶下標(biāo)L, R的符號(hào)Uu UK分別表示U在xs點(diǎn)左邊和右邊的 值。
[0008] 守恒格式的差分解為： 卜（./<0
[0009] U- - -,UU (7 - /); \UR (./> 0
[0010] 其中，下標(biāo)j為網(wǎng)格節(jié)點(diǎn)編號(hào)，仏為\點(diǎn)上的近似U值。
[0011] 為了滿足守恒性，需要下式成立： _2] Uy (X +, ^^) = UL (.τ, ^x,H) + UR (x+ ； ^)-
[0013] 此式和差分格式唯一地確定了XJPUM。因此，數(shù)值間斷必須至少涉及三個(gè)格點(diǎn)（即 兩個(gè)區(qū)間）。目前一些優(yōu)良格式的間斷分辨率已達(dá)到二格至三格。要實(shí)現(xiàn)高分辨率格式的 設(shè)計(jì)，就是在保證格式的相容性、穩(wěn)定性的前提下，在設(shè)計(jì)和構(gòu)造格式時(shí)能夠盡可能減弱其 數(shù)值耗散性，保持?jǐn)?shù)值耗散性優(yōu)勢(shì)。
[0014] 網(wǎng)格生成是數(shù)值模擬的基礎(chǔ)，它實(shí)質(zhì)是物理求解域與計(jì)算求解域的轉(zhuǎn)換。在求解 具有復(fù)雜幾何形狀的流場(chǎng)時(shí)，適當(dāng)?shù)木W(wǎng)格生成是一個(gè)十分重要的問(wèn)題，網(wǎng)格質(zhì)量的好壞直 接影響到計(jì)算結(jié)果的收斂情況及精度。無(wú)結(jié)構(gòu)網(wǎng)格具有復(fù)雜區(qū)域適應(yīng)性好、局部加密靈活 和便于自適應(yīng)的優(yōu)點(diǎn)，能很好地模擬自然邊界及復(fù)雜的水下地形，提高邊界模擬精度。


【發(fā)明內(nèi)容】

[0015] 本發(fā)明所要解決的技術(shù)問(wèn)題是：針對(duì)現(xiàn)有技術(shù)的缺陷，提供一種求解淺水問(wèn)題模 擬間斷水流數(shù)值的方法，，具體是一種求解二維淺水波方程黎曼解的數(shù)值方法，可以應(yīng)用于 模擬潰壩波、涌潮等不可壓間斷水流，進(jìn)而了解間斷水流的水動(dòng)力特征。
[0016] 本發(fā)明為解決上述技術(shù)問(wèn)題采用以下技術(shù)方案：
[0017] 一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，建立模擬間斷水流問(wèn)題的高分辨率 數(shù)值模型，利用所述模型模擬地形和計(jì)算域中帶有大梯度解或間斷的水流現(xiàn)象，具體步驟 如下：
[0018] 步驟一、建立計(jì)算網(wǎng)格，將物理計(jì)算區(qū)域進(jìn)行無(wú)結(jié)構(gòu)三角形網(wǎng)格分割，得到復(fù)數(shù)個(gè) 三角形網(wǎng)格單元，將單一的網(wǎng)格單元作為控制元，其中第i個(gè)三角形控制元記為;
[0019] 步驟二、通過(guò)守恒型物理量反映每個(gè)計(jì)算網(wǎng)格的水動(dòng)力特征，在笛卡爾直角坐標(biāo) 系下，將平面二維淺水運(yùn)動(dòng)守恒形式的方程表示為：
[0020] Ut+WU)x+G^y = S (1);
[0021] 其中，U為守恒量向量；F，G為通量向量；S為源項(xiàng)向量，t為時(shí)間，X和y為笛卡爾 坐標(biāo)系坐標(biāo)，U t，F(xiàn)(U)X，G(U)y分別表示為時(shí)間t和空間X，y方向的偏導(dǎo)數(shù)，具體表示為：
[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表示向量的轉(zhuǎn)置；u、v分別是x、y方向沿水深平均的流速分量，并且為 關(guān)于X，y，t的函數(shù)；h是水深；g是重力加速度；b x、by分別為河床比降在X、y方向上的分 量，b (X，y)為河底高程函數(shù)；
[0027] 求解上述方程的黎曼解，對(duì)方程組（1)，令向量函數(shù)E= (F，G)，V= 為梯 ^ cx cn- J 度算子，將式（1)表示為：
[0028] Ut+ ▽ · E = S (2)；
[0029] 采用步驟一中劃分的三角形網(wǎng)格單元對(duì)計(jì)算區(qū)域進(jìn)行離散，在上對(duì)方程組（2) 進(jìn)行積分，并按逆時(shí)針利用Green公式將面積分化為線積分，得：

【權(quán)利要求】
1. 一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，其特征在于，建立模擬間斷水流問(wèn)題 的高分辨率數(shù)值模型，利用所述模型模擬地形和計(jì)算域中帶有大梯度解或間斷的水流現(xiàn) 象，具體步驟如下： 步驟一、建立計(jì)算網(wǎng)格，將物理計(jì)算區(qū)域進(jìn)行無(wú)結(jié)構(gòu)三角形網(wǎng)格分割，得到復(fù)數(shù)個(gè)三角 形網(wǎng)格單元，將單一的網(wǎng)格單元作為控制元，其中第i個(gè)三角形控制元記為; 步驟二、通過(guò)守恒型物理量反映每個(gè)計(jì)算網(wǎng)格的水動(dòng)力特征，在笛卡爾直角坐標(biāo)系下， 將平面二維淺水運(yùn)動(dòng)守恒形式的方程表示為： Ut+F(U)x+G(U)y = S (1); 其中，U為守恒量向量；F，G為通量向量；S為源項(xiàng)向量，t為時(shí)間，X和y為笛卡爾坐標(biāo) 系坐標(biāo)，Ut，F(xiàn)(U)X，G(U)y分別表示為時(shí)間t和空間X，y方向的偏導(dǎo)數(shù)，具體表示為 ： 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表示向量的轉(zhuǎn)置；u、v分別是x、y方向沿水深平均的流速分量，并且為關(guān)于 X，y，t的函數(shù)；h是水深；g是重力加速度；bx、by分別為河床比降在x、y方向上的分量； 求解上述方程的黎曼解，對(duì)公式（1)，令向量函數(shù)E = (F，G)，▽ 為梯度算 ^ ox oy ) 子，將式⑴表示為： Ut+ ▽ · E = S (2)； 采用步驟一中劃分的三角形網(wǎng)格單元對(duì)計(jì)算區(qū)域進(jìn)行離散，在上對(duì)方程組（2)進(jìn) 行積分，并按逆時(shí)針利用Green公式將面積分化為線積分，得： j 3 UdQ = -Y ? E-nJI+? Sdil (3)； di tfiL·-·· Jn.+ 其中，為時(shí)間方向的常微分算子，cm是面積分微元，di是線積分微元，Lm表示第i di 個(gè)三角形單元Ω i的第m條邊，m = 1，2, 3, nm = (nmx, nmy) = (cos θ，sin θ )為L(zhǎng)m的單位外 法向量，nmx，nmySnm在x、y方向的分量，Θ為外法線方向與水平方向的夾角； 定義單元平均值士由1嚴(yán)，&赫，亂其中，ΙΩ」為網(wǎng)格單元面積； 利用兩點(diǎn)Gauss求積公式求解式（3)中的線積分： I ￡ ·njl?|LmImqEψ nm ⑷； m q=i 其中，I Lm I為三角形邊界線Lm的長(zhǎng)度，G,為三角形單元邊Lm上的Gauss點(diǎn)，(f = u) _ 1 為相應(yīng)的權(quán)系數(shù)，q為咼斯點(diǎn)編號(hào)，?I = 0? ; 構(gòu)造一個(gè)數(shù)值通量格式近似替代E(U(Gq，t))，將三角形網(wǎng)格單元邊界上內(nèi)部和外部 的變量值分別記為隊(duì)、UK，將通過(guò)Gauss點(diǎn)G,的通量記為，其中， UjGq，t)，UjGpt)分別為U在Gauss點(diǎn)三角形網(wǎng)格單元內(nèi)部和外部的值，由此得出控制方 程（3)的有限體積半離散化格式：
步驟三、求解三角形網(wǎng)格單元邊界上的HLL數(shù)值通量# ; 步驟四、采用三階Runge-Kutta時(shí)間離散方法對(duì)時(shí)間方向上的離散進(jìn)行處理，得出三 階精度的時(shí)間離散格式； 步驟五、采用三階WENO格式重構(gòu)方法求解UJG,，t)，UJG,，t); 步驟六、對(duì)每個(gè)網(wǎng)格單元，利用當(dāng)前時(shí)間層上每個(gè)單元的守恒量向量積分平均值，求解 下一時(shí)刻的守恒量向量積分平均值： 經(jīng)過(guò)步驟五對(duì)流體方程時(shí)間方向上的的離散后，再對(duì)流體方程進(jìn)行空間離散，表達(dá)式 為： ^ = L(U)； di 其中，L(U)是空間離散化算子； 設(shè)定第η個(gè)時(shí)間層上的U值為Un，則經(jīng)過(guò)一個(gè)時(shí)間步長(zhǎng)Λ t后，下個(gè)時(shí)間層的近似值 Un+1表示為： U(1) = Un+AtL(Un)；
其中，U(1)、U(2)為中間過(guò)渡向量。
2.如權(quán)利要求1所述的一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，其特征在于，所 述步驟三中，方程的HLL通量表達(dá)式具體如下：
其中，sK分別是波速度的上下確界At = %-a!^，sK = uK+aKqK ; 其中，左右兩邊的波速表示為：? = =^Γ; 設(shè)定qK、f為中間過(guò)渡量，Κ = L，R，具體為：
3.如權(quán)利要求1所述的一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，其特征在于，所 述步驟四包括：針對(duì)三角形網(wǎng)格單元構(gòu)造三角形單元的多項(xiàng)式，具體步驟包括： (501) 線性多項(xiàng)式的構(gòu)造： 基于三角形單元以及相鄰的九個(gè)三角形單元，設(shè)定九個(gè)模板，具體為： 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}； 其中，下標(biāo)i, j, k, ia, ib, ja, jb, ka, kb為不同相鄰單元的編號(hào)； 定義T = U B2 U…U B9為混合模板集，其中U為并集符號(hào)，在每個(gè)小模板上構(gòu)造線 性多項(xiàng)式P(x，y)，使得P(x，y)在小模板所覆蓋的三個(gè)網(wǎng)格單元上的平均值和已知單元Ω。 的平均值相等； (502) 構(gòu)造線性權(quán)： 構(gòu)造如下二次多項(xiàng)式：Q(x，y) = h+t^x+t^y+t^xy+by+bjjy2 ; 其中,131，132，133，134，13 5，136為多項(xiàng)式的待定系數(shù)，在{〇(|，1，」， 1;，^，如0>0# Qka, Qkb}這十個(gè)單元上滿足：, s = 0, i, j, k, ia, ib, ja, jb, ka, kb ;其中 ￡/s 表示單元Ω s上U的平均值， 對(duì)三角形單元Ω ^的每個(gè)Gauss積分點(diǎn)（xG, yG)，設(shè)定線性權(quán)為γ i (1 = 1，. . .，9), γ 1 是由三角形網(wǎng)格的幾何性質(zhì)唯一決定的常數(shù)，權(quán)系數(shù)非負(fù)； 將線性多項(xiàng)式Ρ (X，y)用線性權(quán)Y i合并： 傘，) = Σ#(χ，.ν) (7) i=I y i 滿足下述等式：R(xe, ye) = Q(xe, ye) (8) (503) 構(gòu)造非線性權(quán)和光滑因子： 在線性權(quán)的基礎(chǔ)上構(gòu)造非線性權(quán)，基于一致性，權(quán)系數(shù)和為1，即Σα二】； /~ι 將九個(gè)線性多項(xiàng)式分為四組： /=1 /=1 對(duì)于i邊上的第一個(gè)Gauss點(diǎn)： 第一組含有方程 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,： 得出線性多項(xiàng)式：壞^卜工方巧卜^^為多項(xiàng)式的非負(fù)解； /-I Π 令：) = ^w。+ a2u, + 〇義+α4-- + a5?ift，對(duì)一個(gè)n次多項(xiàng)式p (X，y)，定義光滑
測(cè)量：ψ= Σ ]>Γ'(dXu))2 Ι<--/|<? 其中，a是一個(gè)復(fù)指數(shù)，D是偏導(dǎo)數(shù)，則非線性權(quán)為
其中，-?是第?個(gè)多項(xiàng)式^的光滑度，ε是小正數(shù)。
4. 如權(quán)利要求3所述的一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，其特征在于：所 述ε的取值范圍為（1(Γ6, ΚΓ2)。
5. 如權(quán)利要求4所述的一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，其特征在于：所 述ε的取值為ε = ΚΓ6。
6. 如權(quán)利要求1所述的一種求解淺水問(wèn)題模擬間斷水流數(shù)值的方法，其特征在于：所 述步驟六中，時(shí)間步長(zhǎng)At的具體計(jì)算方式為：
其中，CFL表示CFL條件數(shù)。
【文檔編號(hào)】G06F19/00GK104091065SQ201410315720
【公開(kāi)日】2014年10月8日 申請(qǐng)日期:2014年7月3日 優(yōu)先權(quán)日:2014年7月3日
【發(fā)明者】盧長(zhǎng)娜 申請(qǐng)人:南京信息工程大學(xué)