一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法
【專利摘要】本發(fā)明公開了一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，方法步驟如下，步驟1，構(gòu)建描述坡面流運(yùn)動(dòng)的運(yùn)動(dòng)波方程，并設(shè)置初始和邊界條件；步驟2，在離散速度模型的基礎(chǔ)上構(gòu)建演進(jìn)方程；步驟3，對(duì)Lattice Boltzmann方法的多尺度問題進(jìn)行處理計(jì)算，步驟4，進(jìn)行平衡態(tài)分布函數(shù)的確定；步驟5，對(duì)初始條件和邊界條件進(jìn)行處理；步驟6，通過改進(jìn)的Lattice Boltzmann法求解坡面流運(yùn)動(dòng)方程。與現(xiàn)有技術(shù)相比較，本發(fā)明可應(yīng)用于坡面流運(yùn)動(dòng)方程的求解，通過改進(jìn)多尺度問題的處理方法，有效地提高了坡面流運(yùn)動(dòng)方程的求解精度，解決了城市化進(jìn)程加速的背景下，流域下墊面地形及覆被特性日益復(fù)雜等問題。
【專利說明】
一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及計(jì)算流體動(dòng)力學(xué)領(lǐng)域，尤其涉及一種城市復(fù)雜下墊面條件下的匯流計(jì) 算方法。
【背景技術(shù)】
[0002] 坡面水流是產(chǎn)匯流研究的重要組成部分，在水文學(xué)及應(yīng)用研究領(lǐng)域中占有舉足輕 重的地位；同時(shí)，坡面水流是土壤水蝕過程的主要?jiǎng)恿?，因此?zhǔn)確模擬坡面水流運(yùn)動(dòng)過程， 對(duì)于揭示坡面水流運(yùn)動(dòng)規(guī)律，明確坡面侵蝕機(jī)理，構(gòu)建土壤侵蝕物理過程模擬具有重要意 義。
[0003] 坡面水流運(yùn)動(dòng)屬于流體運(yùn)動(dòng)的范疇。目前對(duì)流體運(yùn)動(dòng)的描述可從宏觀、介觀、微觀 3種觀察尺度進(jìn)行。宏觀連續(xù)介質(zhì)模型將流體假設(shè)為連續(xù)介質(zhì)并用Navier Stokes方程組及 其簡化形式來描述，然后采用某種離散方法進(jìn)行求解。這類方法是目前流體力學(xué)計(jì)算中最 成熟和最成功的方法。但不可否認(rèn)在宏觀層次上描述流體運(yùn)動(dòng)還面臨著挑戰(zhàn)：①對(duì)于一些 復(fù)雜情況，如湍流、多孔介質(zhì)流、跨尺度流動(dòng)等，由于流動(dòng)機(jī)制或邊界條件復(fù)雜，宏觀模型有 很大的局限性。②對(duì)于強(qiáng)非線性的微分方程，數(shù)值穩(wěn)定性是一個(gè)很大的問題。③對(duì)于每一時(shí) 間步長，都需要求解大型線性代數(shù)方程組，對(duì)計(jì)算機(jī)的計(jì)算速度有很高的要求。微觀分子模 型將流體視為由大量分子構(gòu)成的系統(tǒng)，著眼于每個(gè)流體分子的動(dòng)力學(xué)行為，通過對(duì)每個(gè)分 子運(yùn)動(dòng)的刻畫來描述流體的整體運(yùn)動(dòng)情況。由于微觀分子模型基于最基本的分子運(yùn)動(dòng)規(guī) 律，原則上可用于模擬任意的流體系統(tǒng)，但它著眼于每個(gè)流體分子的動(dòng)力學(xué)行為，對(duì)于任一 個(gè)流體系統(tǒng)分子數(shù)目都是非常龐大的，這顯然需要巨大的計(jì)算量和存儲(chǔ)量，在實(shí)際應(yīng)用中 并不可行。介觀動(dòng)理學(xué)模型建立在Boltzmann方程的基礎(chǔ)上，它著眼于流體分子的速度分布 函數(shù)，并認(rèn)為，單個(gè)分子的運(yùn)動(dòng)細(xì)節(jié)并不影響流體的宏觀特性，可通過構(gòu)造一些簡單的演化 規(guī)則，來獲得與物理規(guī)律相符的數(shù)值結(jié)果。介觀動(dòng)理學(xué)模型既可以從微觀角度考察流體分 子的運(yùn)動(dòng)信息，避免了宏觀模型連續(xù)介質(zhì)的假設(shè)，又避免了微觀分子模型分子數(shù)目龐大的 問題。目前，常見的介觀模擬方法有格子氣自動(dòng)機(jī)、Lattice Boltzmann方法、直接蒙特卡羅 方法等。其中，Lattice BGK模型的提出，使得復(fù)雜的碰撞操作轉(zhuǎn)化成一個(gè)簡單的松弛過程， 矩陣由松弛時(shí)間確定，大大簡化了計(jì)算，極大地提高了 La11 i ce Bo 11zmann方法的計(jì)算效 率，在國內(nèi)外引起了廣泛的關(guān)注，但迄今有關(guān)Lattice Boltzmann方法應(yīng)用于坡面水流運(yùn)動(dòng) 的研究鮮見報(bào)道。本發(fā)明運(yùn)用Lattice Boltzmann方法，求解坡面流運(yùn)動(dòng)方程，并采用人工 模擬降雨坡面水流實(shí)驗(yàn)驗(yàn)證其有效性。

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

[0004] 本發(fā)明的目的就在于提供一種解決上述問題，基于Lattice Boltzmann提高城市 復(fù)雜下墊面條件下流域產(chǎn)匯流計(jì)算的精度以及效率的城市復(fù)雜下墊面條件下的匯流計(jì)算 方法。
[0005] 為了實(shí)現(xiàn)上述目的，本發(fā)明采用的技術(shù)方案是:一種城市復(fù)雜下墊面條件下的匯 流計(jì)算方法，方法步驟如下，
[0006] 步驟1，構(gòu)建描述坡面流運(yùn)動(dòng)的運(yùn)動(dòng)波方程，并設(shè)置初始和邊界條件；
[0007] 步驟2,在離散速度模型的基礎(chǔ)上構(gòu)建Lattice Boltzmann方法的演進(jìn)方程；
[0008] 步驟3,改進(jìn)Lattice Boltzmann方法多尺度問題的處理方法，具體為：
[0009] 傳統(tǒng)Lattice Boltzmann方法在進(jìn)行多尺度處理時(shí)，通常是引入兩個(gè)時(shí)間尺度和 一個(gè)空間尺度，但是這種模型對(duì)邊界格式非常敏感，分布函數(shù)經(jīng)常出現(xiàn)負(fù)值，違背了物理規(guī) 律，難以得到正確結(jié)果，但如果對(duì)空間坐標(biāo)不進(jìn)行多重尺度化處理，只對(duì)時(shí)間作多重尺度化 處理則可以解決此問題。
[0010] 對(duì)D1Q3模型，如空間坐標(biāo)不進(jìn)行多尺度處理，而對(duì)時(shí)間坐標(biāo)引入3個(gè)時(shí)間尺度：tk = ekt(k = 0、l、2)，則可得到時(shí)間和空間坐標(biāo)的多尺度表達(dá)的微分形式為
[0011: (28)
[0012]式中ε為Knudsen數(shù)，定義為分子平均自由程和流動(dòng)的宏觀特征長度之比，是衡量 流體稀薄程度的無因次量。
[0013] 將演進(jìn)方程的fa(x+ea Δ t，ea，t+ Δ t)作二階Taylor展開，并用Knudsen數(shù)ε代替時(shí) 間步長Δ t,整理分析可得剞3個(gè)不同時(shí)間R庶的富散Roltzmann方程，即
[0014]
[0015] ,、 (30)
[0016] (31)
[0017] 步驟4，平衡態(tài)分布函數(shù)的確定，具體為：
[0018] 從數(shù)學(xué)角度出發(fā)，以多尺度分析為手段通過待定系數(shù)法來確定平衡態(tài)分布函數(shù)， 離散的粒子速度分布函數(shù)采用Chapmann Enskog展開為
[0019] fa = fa(0)+efa(1)+e2fa(2)+O(e 3) (32)
[0025] 對(duì)式(3)各項(xiàng)取零階矩，得
[0020]由于在Chapmann Enskog展開中假設(shè)局部范圍內(nèi)已接近平衡態(tài)，故對(duì)式(8)兩邊求 和，有
[0021]
[0022]
[0023]
[0024]
[0026](36)
[0027]坡面流運(yùn)動(dòng)的基本方程都是修正的Burgers方程的特例。修正的Burgers方程是一 個(gè)非線性對(duì)流-擴(kuò)散方程，形式為[0028](37)[0029] 式中F為源項(xiàng)，即外力作用項(xiàng);u為變量，代表水深、濃度、含水率、流量等;β為阻力 參數(shù);λ、1為參數(shù);m反映流體流動(dòng)狀態(tài)的參數(shù)。[0030] 將忒（12)與忒講桿相比，容易得到[0031] (38)[0032] (39)[0033] 將式（14)和式（15)代入式（12)和式（13)中，可得[0034][0035][0036][0037][0038] 這樣，通過求取各階矩的方式就建立了宏觀量與平衡態(tài)分布函數(shù)之間的聯(lián)系。為 了還原修正的Burgers方程，需對(duì)不同時(shí)間尺度的離散Boltzmann方程在不同方向上進(jìn)行多 尺度粘合(即將時(shí)間尺度t0、tl、t2還原回時(shí)間尺度t)。[0039] 式(9)兩邊同時(shí)乘以ε，并與式(3)兩邊相加，然后對(duì)3個(gè)方向求和，可得[0040][0041][0042][0043]比較式(20)和式（13)可看出[0044] F = e(g〇+gi+g2) (45)[0045] 由式(21)可知，無論ga取什么值，只要ga3個(gè)方向取值之和與ε的乘積等于F即可，本 文取平均值，即[0046] ga = F/(3e)(a = 〇,U2) (46)[0047] 另外，式(20)右:'多了_ε和增加了〇(ε2)項(xiàng)，其中〇(ε2)的大小由弛 豫時(shí)間τ來控制，而通過宏； S新選取可實(shí)現(xiàn)_ε的消除。令




[_]
(47)
[0049] 式中，k為可調(diào)參數(shù)。
[0050]當(dāng)預(yù)先設(shè)定弛豫時(shí)間τ和Knudsen數(shù)后，k的計(jì)算公式為 [0051]
(48)
[0052] 經(jīng)過這樣的處理后，即可還原出修正的Burgers方程。對(duì)于D1Q3模型，3個(gè)方向的平 衡態(tài)分布函數(shù)可由式（14)、（15)和（18)確定。由于el = -c、e0 = 0、e2 = c，故D1Q3模型對(duì)應(yīng)的 修正的Burgers方程的平衡態(tài)分布函數(shù)為
[0053] (49)
[0054]步驟5，初始條件和邊界條件的處理，具體為：
[0057]式中f(n)為經(jīng)過迀移后的粒子速度分布函數(shù)。[0058]下邊界結(jié)點(diǎn)(M = N)上的分布函數(shù)為
[0055] 上邊界結(jié)點(diǎn)上的分布函數(shù)用平衡態(tài)分布函數(shù)代替，以宏觀上邊界條件作為限制條 件;通過分布函數(shù)外推格式來確定下邊界結(jié)點(diǎn)上的分布函數(shù)。則對(duì)于D1Q3模型，上邊界結(jié)點(diǎn) (M=O)上的分布函數(shù)為
[0056] (50)
[0059] (51)
[0060] 式中f(rrf)為粒子迀移前的速度分布函數(shù)，f(N-l)、f(N_2)為第N-1、N_2個(gè)節(jié)點(diǎn)上的 分布函數(shù)。
[0061] 步驟6,改進(jìn)的Lattice Boltzmann法求解坡面流運(yùn)動(dòng)方程的具體步驟如下：
[0062] 首先把空間劃分成均勻的網(wǎng)格(格子），選擇合適的弛豫時(shí)間τ，然后根據(jù)初始時(shí)刻 的水深和速度，求得每個(gè)結(jié)點(diǎn)各方向的平衡態(tài)分布函數(shù)f a，初始時(shí)刻各結(jié)點(diǎn)的分布函數(shù)用 平衡態(tài)分布函數(shù)代替，即式(25);使Latti ce Bo11zmann離散方程(7)在網(wǎng)格結(jié)點(diǎn)的各方向 上進(jìn)行演進(jìn)(迀移和碰撞過程），即令fa按式(7)在網(wǎng)格上迭代求解。對(duì)分布函數(shù)求各階矩， 計(jì)算出該時(shí)刻所對(duì)應(yīng)的宏觀量水深等。具體地來說，碰撞過程采用tn時(shí)刻的水深h和速度V 計(jì)算各結(jié)點(diǎn)各方向的平衡分布函數(shù)量/:，然后根據(jù)式(7)計(jì)算tn時(shí)刻碰撞后各結(jié)點(diǎn)各方向 的分布函數(shù)fa;迀移過程是各結(jié)點(diǎn)各方向的分布函數(shù)fV沿對(duì)應(yīng)的方向迀移到新的結(jié)點(diǎn)，得到 tn+Ι時(shí)刻的分布函數(shù)fa。更新水深h和速度V;再進(jìn)行tn+Ι時(shí)刻后碰撞和迀移過程的演算，跳 回循環(huán)計(jì)算，直到滿足誤差為止。
[0063] 作為優(yōu)選，所述步驟1中構(gòu)建的坡面流運(yùn)動(dòng)波方程為：
[0064]
(52)
[0065] 具中，So = Sin^，r = p_i。
[0066] 所描述步驟1中初始條件和邊界條件設(shè)置為
[0067]
(53)
[0068]式中q為單款流量，m2/s ;h為水深m; X為坡面某點(diǎn)至坡頂?shù)木嚯x，m; t為時(shí)間，s ;n為 曼寧糙率系數(shù);SO為坡面坡度；Θ為坡面傾角;r為凈雨強(qiáng)度，p為降雨強(qiáng)度，m/s;i為下滲率， m/s〇
[0069]作為優(yōu)選，所述步驟2是基于D1Q3離散速度模型建立Lattice Boltzmann方法的演 進(jìn)方程。D1Q3表示空間維數(shù)取1，離散速度方向數(shù)取S13Lattice Boltzmann方法的演進(jìn)方程 為
[0070]
：54)
[0071]式中，α為離散的粒子運(yùn)動(dòng)的方向序號(hào);ea為離散的粒子在第α方向上的運(yùn)動(dòng)速度， m/s; △ t為時(shí)間步長，s;x為粒子的位置矢量;fa(x，ea，t)為尚散的粒子在第a方向上的速度 分布函數(shù)；Qa(f( x，ea，t))為第α方向上離散的粒子的碰撞項(xiàng)；τ為馳豫時(shí)間（〇 < τ )，s ; /，(^^…為平衡態(tài)分布函數(shù)"^^⑴&⑴為源項(xiàng)…為粒子運(yùn)動(dòng)方向的總數(shù)志離散速 度模型有關(guān)。
[0072]與現(xiàn)有技術(shù)相比，本發(fā)明的優(yōu)點(diǎn)在于:本發(fā)明基于Lattice Boltzmann方法，將其 應(yīng)用于坡面流運(yùn)動(dòng)方程的求解，通過改進(jìn)多尺度問題的處理方法，有效地提高了坡面流運(yùn) 動(dòng)方程的求解精度。有效解決了城市化進(jìn)程加速的背景下，流域下墊面地形及覆被特性日 益復(fù)雜，傳統(tǒng)坡面匯流計(jì)算方法在計(jì)算精度和計(jì)算時(shí)效上已經(jīng)很難滿足當(dāng)前城市突發(fā)暴雨 洪水計(jì)算需求的問題。在流域產(chǎn)匯流模型構(gòu)建及其求解、土壤侵蝕物理過程模擬等方面將 具有廣闊的前景。
【附圖說明】
[0073]圖1為本發(fā)明的計(jì)算流程圖。
[0074]圖2為實(shí)測(cè)坡面徑流與計(jì)算坡面徑流的比較。
[0075] 圖3為格子Boltzmann方法和Prei ssmann4點(diǎn)隱式格式法在坡面流運(yùn)動(dòng)過程中應(yīng)用 的對(duì)比分析。
[0076] 圖4為格子Boltzmann方法和Prei ssmann4點(diǎn)隱式格式法在坡面流運(yùn)動(dòng)過程中應(yīng)用 的計(jì)算相對(duì)誤差。
【具體實(shí)施方式】
[0077]下面將對(duì)本發(fā)明作進(jìn)一步說明。
[0078]實(shí)施例:在西北農(nóng)林科技大學(xué)黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點(diǎn)實(shí)驗(yàn)室構(gòu)建 3個(gè)可變坡實(shí)驗(yàn)槽，每個(gè)實(shí)驗(yàn)槽長、寬和高分為2m、0.55m和0.3m。實(shí)驗(yàn)槽中，設(shè)有地面徑流出 口，壤中流出口和地下徑流3個(gè)出口。3個(gè)徑流出口位置分別設(shè)在距槽底0.35、0.175和0cm。 實(shí)驗(yàn)槽的下端配置有1個(gè)"V"型集流裝置，用于定期采集徑流樣。
[0079] 土樣為重壤土，取自陜西楊凌的實(shí)驗(yàn)田，采用分層取土法。分層填土完畢后，用塑 料薄膜覆蓋土槽表面，以防異物落入。3個(gè)實(shí)驗(yàn)槽的不同之處在于下墊面情況不同。其中實(shí) 驗(yàn)槽1里面種植有紫花苜蓿草，實(shí)驗(yàn)槽2里面插有人工大孔隙，實(shí)驗(yàn)槽3為裸土坡面。
[0080]人工降雨系統(tǒng)采用黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點(diǎn)實(shí)驗(yàn)室的側(cè)噴區(qū)自動(dòng) 模擬降雨系統(tǒng)，該降雨系統(tǒng)實(shí)際降雨高度為16m，雨滴達(dá)到的終點(diǎn)速度滿足天然降雨特性， 降雨均勻度高達(dá)85 %以上，降雨強(qiáng)度可在15~260mm/h范圍內(nèi)分級(jí)可調(diào)，可持續(xù)降雨 lOOmin。將3個(gè)土槽的坡度設(shè)定在5°，共進(jìn)行了3場(chǎng)次降雨，降雨強(qiáng)度為30、60、120mm/h。
[0081 ] 產(chǎn)流計(jì)算采用VIMAC模型，坡面匯流計(jì)算方法則采用Lattice Boltzmann法，其結(jié) 果見圖2所示。圖中Q為坡面徑流量，p為降雨強(qiáng)度。
[0082]對(duì)照?qǐng)D2，2 0 0 71013場(chǎng)次槽1、槽2和槽3實(shí)驗(yàn)資料的地面徑流深誤差分別為：-1.4721 %、-7·3021 %和1.0156%，確定性系數(shù)分別為：0.57、0·51 和0.68;20071027場(chǎng)次槽3 的地面徑流深誤差分別為：-1.2596%、-10· 6626%、2· 3351 %，確定性系數(shù)分別為：0· 88、 0.80和0.93; 20071102場(chǎng)次槽3地面徑流深誤差分別為：1.5068%、4.2286%、1.5149%，確 定性系數(shù)分別為:〇.95、0.79和0.91。即3場(chǎng)次實(shí)驗(yàn)資料的地面徑流深誤差均在±11%以內(nèi)， 確定性系數(shù)均大于0.51，其中大于0.70的占67 %，最大的達(dá)0.95。這說明該Lattice Bo 11 zmann方法可成功應(yīng)用于坡面流運(yùn)動(dòng)方程的求解。
[0083]另外，通過理想算例，以解析解為標(biāo)準(zhǔn)，對(duì)格子Boltzmann方法與應(yīng)用廣泛的 Preissmann 4點(diǎn)隱式差分法的計(jì)算精度進(jìn)行比較分析。
[0084] 理想化算例構(gòu)建為:假定一規(guī)則矩形不透水坡面，坡面長L = 50m，糙率n = 0.02,坡 度5 = 0.01，研究在降雨強(qiáng)度為3〇1111]1/11、降雨歷時(shí)為1〇1]1；[11、退水階段歷時(shí)401]1；[11情況下的坡 腳水深和流量隨時(shí)間的變化以及坡面水深和流量沿坡長變化情況。
[0085]對(duì)Preissmann 4點(diǎn)隱式差分法，根據(jù)其穩(wěn)定性條件，取加權(quán)系數(shù)Θ = 〇 . 55，空間步 長Δ X = Im，時(shí)間步長Δ t = Is;采用格子BoItzmann方法計(jì)算坡面水流時(shí)，采用相同空間步 長和時(shí)間步長，即空間步長Δ X = lm，時(shí)間步長Δ t = Is，τ是反映粒子分布函數(shù)趨于平衡態(tài) 的時(shí)間（0<τ<2)，這里取τ = 1.2;因?yàn)槠旅娌煌杆越涤険p失可以忽略。
[0086] 對(duì)照?qǐng)D3,由圖3a和3b可以看出，格子Boltzmann方法的D1Q5速度模型和 Preissmann4點(diǎn)隱式差分法與解析解的計(jì)算結(jié)果幾乎一致。由圖3c~3d可以看出，坡頂水深 始終是零，隨著離坡頂距離增加，水深逐漸變厚，并在后端形成一平頂。但由于在坡腳水流 也在不斷流出，水流厚度并不會(huì)無限增大，在一段時(shí)間后達(dá)到平衡，平頂消失，水深自上而 下變厚，但其分布趨于穩(wěn)定，見圖3e和3f。格子Bol tzmann方法的D1Q5速度模型計(jì)算結(jié)果表 明，在t = 7.OSmin時(shí)，平頂區(qū)域完全消失，不再隨時(shí)間而變。這意味著達(dá)到動(dòng)態(tài)平衡時(shí)間點(diǎn)， 從外界輸入的水量和從坡腳流走的水量相等。在以后的時(shí)間里，只要外界輸入項(xiàng)不變，坡面 上水深分布形狀就保持不變，坡腳流量也為一常數(shù)。所以t = 5min時(shí)，單位時(shí)間從外界降雨 輸入的水量大于從坡腳流走的水量，即圖3c和3d刻畫的是未達(dá)到平衡時(shí)間之前的坡面水深 與單寬流量沿坡長變化;t= IOmin時(shí)，單位時(shí)間從外界降雨輸入的水量等于從坡腳流走的 水量，即圖3e和3f刻畫的是達(dá)到平衡時(shí)間之后的坡面水深與單寬流量沿坡長變化。從圖3c 和3d可以看出，PreiSSmann4點(diǎn)隱式差分法在平衡位置以后有稍許數(shù)值振蕩，而格子 Boltzmann方法的D1Q5速度模型自始至終未出現(xiàn)數(shù)值振蕩；從圖3e和圖3f可以看出，格子 13〇]^211^1111方法的01( >)5速度模型和?161881]^11114點(diǎn)隱式差分法與解析解結(jié)果非常一致。 [0087] 為了對(duì)格子Boltzmann方法的D1Q5速度模型和Preissmann4點(diǎn)隱式差分法的計(jì)算 精度作進(jìn)一步分析，將解析解作為標(biāo)準(zhǔn)，采用相對(duì)誤差的概念，確定格子Bol tzmann方法的 Dl Q5模型的相對(duì)誤差和Pre i ssmann4點(diǎn)隱式差分法的相對(duì)誤差，結(jié)果如圖4所示。
[0088]對(duì)照?qǐng)D4,由圖4a和4b可以看出，對(duì)于模擬坡腳斷面水深和單寬流量過程，格子 Bol tzmann方法的D1Q5模型的計(jì)算精度好于Prei ssmann4點(diǎn)隱式差分法，尤其在退水階段。 具體來說，對(duì)于坡腳水深，格子Bol tzmann方法的D1Q5模型的相對(duì)誤差在0 %~2.89 %范圍 之間，最大相對(duì)誤差點(diǎn)出現(xiàn)在t = 24min時(shí);Preissmann4點(diǎn)隱式差分法的相對(duì)誤差在0%~ 41.08%范圍之間，在七=25111；[11之前，相對(duì)誤差控制在0%~1.46%之間，在七=251]1；[11之后相 對(duì)誤差呈現(xiàn)指數(shù)增加的趨勢(shì)。對(duì)于坡腳單寬流量，格子Boltzmann方法的D1Q5模型的相對(duì)誤 差在0 %~4.77%范圍之間；Prei ssmann4點(diǎn)隱式差分法的相對(duì)誤差在0 %~77.46 %范圍之 間，其相對(duì)誤差的變化趨勢(shì)與坡腳水深的相對(duì)誤差變化趨勢(shì)相似。由圖4c到4f可以看出，對(duì) 于坡面水深和坡面單寬流量，兩種方法的相對(duì)誤差均在5%以內(nèi)；對(duì)于達(dá)到平衡時(shí)間之前的 坡面水深和坡面單寬流量，PreissmanM點(diǎn)隱式差分法在平衡位置以后出現(xiàn)誤差振蕩現(xiàn)象， 格子Bo I tzmann方法的Dl Q5速度模型在X = Im點(diǎn)上出現(xiàn)了較大的相對(duì)誤差，其坡面水深相對(duì) 誤差達(dá)到2.83 %，坡面單寬流量相對(duì)誤差達(dá)到4.77 %，在X = 28m點(diǎn)上出現(xiàn)另一個(gè)相對(duì)誤差 峰值，其坡面水深相對(duì)誤差達(dá)到1.88%，坡面單寬流量相對(duì)誤差達(dá)到3.12%，但總的來說， 格子Boltzmann方法的D1Q5速度模型的坡面水深和坡面單寬流量的平均相對(duì)誤差（分別為 0.30%和0.49%)小于PreissmanM點(diǎn)隱式差分法的坡面水深和坡面單寬流量的平均相對(duì) 誤差(分別為0.33%和0.56%);對(duì)于達(dá)到平衡時(shí)間之后的坡面水深和坡面單寬流量，格子 Boltzmann方法的D1Q5速度模型的坡面水深和坡面單寬流量的平均相對(duì)誤差（分別為 0.26%和0.38%)大于PreissmanM點(diǎn)隱式差分法的坡面水深和坡面單寬流量的平均相對(duì) 誤差(分別為0.05%和0.08%)，其主要原因是格子Boltzmann方法的D1Q5速度模型在X= Im 點(diǎn)上出現(xiàn)了較大的相對(duì)誤差（其坡面水深相對(duì)誤差達(dá)到2.83%，坡面單寬流量相對(duì)誤差達(dá) 到4.76%) 0
[0089]由于坡面匯流計(jì)算是流域產(chǎn)匯流計(jì)算的基礎(chǔ)，因此可以預(yù)計(jì)Lattice Boltzmann 方法在流域產(chǎn)匯流模型構(gòu)建及其求解中具有廣闊的前景。
[0090] 本發(fā)明首次將Lattice Boltzmann方法應(yīng)用于坡面水流運(yùn)動(dòng)的研究，并改進(jìn)多尺 度問題的處理方法，有效地提高了坡面流運(yùn)動(dòng)方程的求解精度。傳統(tǒng)Lattice Boltzmann方 法在進(jìn)行多尺度處理時(shí)，通常是引入兩個(gè)時(shí)間尺度和一個(gè)空間尺度，但是這種模型對(duì)邊界 格式非常敏感，分布函數(shù)經(jīng)常出現(xiàn)負(fù)值，違背了物理規(guī)律，難以得到正確結(jié)果。本發(fā)明對(duì)空 間坐標(biāo)不進(jìn)行多重尺度化處理，只對(duì)時(shí)間作多重尺度化處理，巧妙地解決了此問題。通過變 坡土槽試驗(yàn)檢驗(yàn)，本發(fā)明可成功應(yīng)用于坡面流運(yùn)動(dòng)方程的求解，計(jì)算所得的各場(chǎng)次降雨的 地面徑流深誤差均在± 11 %以內(nèi)。通過理想算例，將改進(jìn)的Lattice Boltzmann方法與應(yīng)用 廣泛的Preissmann四點(diǎn)隱式差分法進(jìn)行對(duì)比。結(jié)果表明，對(duì)于模擬坡腳斷面水深和單寬流 量過程以及達(dá)到平衡時(shí)間之前的坡面水深和坡面單寬流量，改進(jìn)的Lattice Boltzmann方 法的計(jì)算精度均高于Pre i ssmann四點(diǎn)隱式差分法。
[0091] 坡面匯流計(jì)算是流域產(chǎn)匯流計(jì)算的基礎(chǔ)，隨著城市化進(jìn)程的加速，流域下墊面地 形及覆被特性日益復(fù)雜，傳統(tǒng)坡面匯流計(jì)算方法在計(jì)算精度和計(jì)算時(shí)效上已經(jīng)很難滿足當(dāng) 前城市突發(fā)暴雨洪水計(jì)算需求?？梢灶A(yù)計(jì)本發(fā)明改進(jìn)的Lattice Boltzmann方法在流域產(chǎn) 匯流模型構(gòu)建及其求解、土壤侵蝕物理過程模擬等方面將具有廣闊的前景。
[0092] 以上對(duì)本發(fā)明所提供的一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法進(jìn)行了詳盡 介紹，本文中應(yīng)用了具體個(gè)例對(duì)本發(fā)明的原理及實(shí)施方式進(jìn)行了闡述，以上實(shí)施例的說明 只是用于幫助理解本發(fā)明的方法及其核心思想；同時(shí)，對(duì)于本領(lǐng)域的一般技術(shù)人員，依據(jù)本 發(fā)明的思想，在【具體實(shí)施方式】及應(yīng)用范圍上均會(huì)有改變之處，對(duì)本發(fā)明的變更和改進(jìn)將是 可能的，而不會(huì)超出附加權(quán)利要求所規(guī)定的構(gòu)思和范圍，綜上所述，本說明書內(nèi)容不應(yīng)理解 為對(duì)本發(fā)明的限制。
【主權(quán)項(xiàng)】
1. 一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，其特征在于:方法步驟如下， 步驟1，構(gòu)建描述坡面流運(yùn)動(dòng)的運(yùn)動(dòng)波方程，并設(shè)置初始和邊界條件； 步驟2，在離散速度模型的基礎(chǔ)上構(gòu)建Lattice Bo 1 tzmann方法的演進(jìn)方程； 步驟3,對(duì)Lattice Boltzmann方法的多尺度問題進(jìn)行處理計(jì)算， 對(duì)D1Q3模型，空間坐標(biāo)不進(jìn)行多尺度處理，而對(duì)時(shí)間坐標(biāo)引入3個(gè)時(shí)間尺度：tk= εΗ (k =〇、1、2)，則可得到時(shí)間和空間坐標(biāo)的多尺度表達(dá)的微分形式為(1) 式中（1) ε為Knudsen數(shù)，定義為分子平均自由程和流動(dòng)的宏觀特征長度之比，是衡量流 體稀薄程度的無因次量； 將演進(jìn)方程的fa(x+ea Δ t，ea，t+ Δ t)作二階化ylor展開，并用Knudsen數(shù)e代替時(shí)間步 長At,整理分析可得到3個(gè)不同時(shí)間尺度的離散Boltzmann方程，即步驟4,進(jìn)行平衡態(tài)分布函數(shù)的確定，從數(shù)學(xué)角度出發(fā)，W多尺度分析為手段通過待定 系數(shù)法來確定平衡態(tài)分布函數(shù)； 步驟5,對(duì)初始條件和邊界條件進(jìn)行處理，上邊界結(jié)點(diǎn)上的分布函數(shù)用平衡態(tài)分布函數(shù) 代替，W宏觀上邊界條件作為限制條件;通過分布函數(shù)外推格式來確定下邊界結(jié)點(diǎn)上的分 布函數(shù)； 步驟6，通過改進(jìn)的Lattice Boltzmann法求解坡面流運(yùn)動(dòng)方程。2. 根據(jù)權(quán)利要求1所述的一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，其特征在于:所 述步驟1中，構(gòu)建的坡面流運(yùn)動(dòng)波方程為(5) 其中，S〇 = sin 目，r = p-i; 所描述步驟1中初始條件和邊界條件設(shè)置為(6) 式中q為單款流量，m2/s ;h為水深m;x為坡面某點(diǎn)至坡頂?shù)木嚯x，m; t為時(shí)間，S ;n為曼寧 糖率系數(shù);SO為坡面坡度;Θ為坡面傾角;r為凈雨強(qiáng)度，P為降雨強(qiáng)度，m/s; i為下滲率，m/s。3. 根據(jù)權(quán)利要求1所述的一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，其特征在于:所 述步驟2是基于D1Q3離散速度模型建立Lattice Boltzmann方法的演進(jìn)方程;D1Q3表示空間 維數(shù)取1，離散速度方向數(shù)取3;Lattice Boltzmann方法的演進(jìn)方程為式中，α為離散的粒子運(yùn)動(dòng)的方向序號(hào)；ea為離散的粒子在第α方向上的運(yùn)動(dòng)速度，m/s; At為時(shí)間步長，s;x為粒子的位置矢量;fa(x，ea，t)為離散的粒子在第α方向上的速度分布 函數(shù)；〇。作^，6。，*))為第〇方向上離散的粒子的碰撞項(xiàng)；1為馳豫時(shí)間（〇《1《〇,3; /。6心，6。，《)為平衡態(tài)分布函數(shù)^〇苗。^，6。山為源項(xiàng);6為粒子運(yùn)動(dòng)方向的總數(shù)，與離散速 度板型有關(guān)。4.根據(jù)權(quán)利要求1所述的一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，其特征在于:步 驟4中，離散的粒子速度分布函數(shù)采用化apmann化skog展開為(8) 由于在化apmann Enskog展開中假設(shè)局部范圍內(nèi)已接近平衡態(tài)，故對(duì)式(8)兩邊求和， 有坡面流運(yùn)動(dòng)的基本方程都是修正的Burgers方程的特例，修正的Burgers方程是一個(gè)非 線性對(duì)流-擴(kuò)散方程，形式為(13) 式中F為源項(xiàng)，即外力作用項(xiàng);U為變量，代表水深、濃度、含水率、流量等;β為阻力參數(shù)； 入、1為參數(shù);m反映流體流動(dòng)狀態(tài)的參數(shù)； 將式(12)與式(13)進(jìn)行相比，容易得到將式(14)和式(15)代入式(12)和式(13)中，可得通過求取各階矩的方式建立了宏觀量與平衡態(tài)分布函數(shù)之間的聯(lián)系，為了還原修正的 Burgers方程，對(duì)不同時(shí)間尺度的離散Boltzmann方程在不同方向上進(jìn)行多尺度粘合，即將 時(shí)間尺度t0、tl、t2還原回時(shí)間尺度t; 式(9)兩邊同時(shí)乘Κε，并與式(3)兩邊相加，然后對(duì)3個(gè)方向求和，可得比較式(20)和式（13)可看出 F=e(gG+gi+g2) (21) 由式(21)可知，無論ga取什么值，只要ga3個(gè)方向取值之和與ε的乘積等于F即可，本文取 平均值，即 邑α二F/(3ε)(日=0、1、2) (22) 式(20)中0(ε2)的大小由弛豫時(shí)間τ來控制，而通過宏觀量均重新選取可實(shí)現(xiàn)-ε 的消除;令(23) 式中，k為可調(diào)參數(shù)； 當(dāng)預(yù)先設(shè)定弛豫時(shí)間τ和Knudsen數(shù)后，k的計(jì)算公式為(24) 通過上述處理，還原出修正的Burgers方程;對(duì)于D1Q3模型，3個(gè)方向的平衡態(tài)分布函數(shù) 可由式（14)、（15)和（18)確定；由于61 = -。、6〇 = 0、62 = (3，故0193模型對(duì)應(yīng)的修正的 Burgers方程的平衡態(tài)分布函數(shù)為(25) ο5. 根據(jù)權(quán)利要求4所述的一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，其特征在于:步 驟5中，對(duì)于D1Q3模型，上邊界結(jié)點(diǎn)(M=0)上的分布函數(shù)為(26) 式中fW為經(jīng)過遷移后的粒子速度分布函數(shù)； 下邊界結(jié)點(diǎn)(M=N)上的分布函數(shù)為(27) 式中為粒子遷移前的速度分布函數(shù)，f(N-l)、f(N-2)為第N-UN-2個(gè)節(jié)點(diǎn)上的分布 函數(shù)。6. 根據(jù)權(quán)利要求5所述的一種城市復(fù)雜下墊面條件下的匯流計(jì)算方法，其特征在于:步 驟6中，求解坡面流運(yùn)動(dòng)方程具體步驟如下，首先把空間劃分成均勻的網(wǎng)格，選擇弛豫時(shí)間 τ，然后根據(jù)初始時(shí)刻的水深和速度，求得每個(gè)結(jié)點(diǎn)各方向的平衡態(tài)分布函數(shù)fa,初始時(shí)刻 各結(jié)點(diǎn)的分布函數(shù)用平衡態(tài)分布函數(shù)代替，即式(25);使Lattice Boltzmann離散方程(7) 在網(wǎng)格結(jié)點(diǎn)的各方向上進(jìn)行演進(jìn)(遷移和碰撞過程），即令fa按式(7)在網(wǎng)格上迭代求解;對(duì) 分布函數(shù)求各階矩，計(jì)算出該時(shí)刻所對(duì)應(yīng)的宏觀量水深等;具體地來說，碰撞過程采用tn時(shí) 刻的水深h和速度V計(jì)算各結(jié)點(diǎn)各方向的平衡分布函數(shù)量/。6"，然后根據(jù)式(7)計(jì)算tn時(shí)刻碰 撞后各結(jié)點(diǎn)各方向的分布函數(shù)fa;遷移過程是各結(jié)點(diǎn)各方向的分布函數(shù)fa沿對(duì)應(yīng)的方向遷 移到新的結(jié)點(diǎn)，得到tn+1時(shí)刻的分布函數(shù)fa;更新水深h和速度V;再進(jìn)行tn+1時(shí)刻后碰撞和 遷移過程的演算，跳回循環(huán)計(jì)算，直到滿足誤差為止。
【文檔編號(hào)】G06F19/00GK105844076SQ201610131568
【公開日】2016年8月10日
【申請(qǐng)日】2016年3月9日
【發(fā)明人】馮杰, 芮孝芳, 楊志勇, 楊濤, 劉寧寧, 何祺勝, 于贏東, 師鵬飛, 王開, 王興勇, 王鵬, 王超
【申請(qǐng)人】中國水利水電科學(xué)研究院, 河海大學(xué)