平面二維水流數(shù)學(xué)模型的建立與應(yīng)用

平面二維水流數(shù)學(xué)模型的建立與應(yīng)用

在計(jì)算二維水流運(yùn)動(dòng)的數(shù)值時(shí)，很難解決區(qū)域邊界的復(fù)雜性。如何處理不規(guī)則邊界，直接影響模型計(jì)算的可靠性和精度。傳統(tǒng)的笛卡爾測(cè)量方法具有簡(jiǎn)單的網(wǎng)格校正原理，適合計(jì)算網(wǎng)格位移，但由于限制因素易于形成鋸齒限制，邊界值調(diào)整的精度不高。為了提高該方法對(duì)非規(guī)則區(qū)域的適應(yīng)性，本工作介紹了一種新的不規(guī)則邊界處理方法。該方法不僅保留了傳統(tǒng)笛卡爾坐標(biāo)方法的簡(jiǎn)單、正交、容易形成計(jì)算網(wǎng)格和計(jì)算效率的特點(diǎn)，而且改進(jìn)了傳統(tǒng)鋸齒限制精度差的缺點(diǎn)。它靠近實(shí)際情況，使邊境地區(qū)的洪水分布變得合理。在復(fù)雜邊界分洪區(qū)的洪水計(jì)算中，具有很高的精度。斜對(duì)角笛卡爾方法基本原則就是按就近逼近原則,把距實(shí)際水邊線(xiàn)較近的網(wǎng)格節(jié)點(diǎn)連成線(xiàn)以擬合實(shí)際水邊線(xiàn).生成網(wǎng)格時(shí),根據(jù)邊界是否固定,可采取不同的判定方法.當(dāng)邊界為固定邊界時(shí),可把固定邊界離散為若干離散節(jié)點(diǎn),確定各離散節(jié)點(diǎn)坐標(biāo),按距離最小原則選定固定邊界上各離散節(jié)點(diǎn)的逼近網(wǎng)格節(jié)點(diǎn),然后依序連接各網(wǎng)格節(jié)點(diǎn)即形成擬合邊界;對(duì)于動(dòng)邊界,由于邊界線(xiàn)變化不定,則可借助網(wǎng)格節(jié)點(diǎn)間水深的干濕變化和水深值,參照相應(yīng)網(wǎng)格節(jié)點(diǎn)間的地形高差,判定距離實(shí)際水邊線(xiàn)較近的網(wǎng)格節(jié)點(diǎn),合并縱橫向的擬合網(wǎng)格節(jié)點(diǎn),形成擬合瞬時(shí)水邊線(xiàn).根據(jù)斜對(duì)角笛卡爾方法基本原則,即以距交點(diǎn)最近為原則,選取對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)作為實(shí)際水邊線(xiàn)的逼近點(diǎn).1模型介紹1.1單寬流量n平面二維水流數(shù)學(xué)模型的基本方程組由水流連續(xù)性方程和水流運(yùn)動(dòng)方程組成.水流連續(xù)性方程:?Z?t+?M?x+?N?y=0(1)?Ζ?t+?Μ?x+?Ν?y=0(1)水流運(yùn)動(dòng)方程:?M?t+?uM?x+?vM?y=?gh?Z?x?gn2uu2+v2√h1/3(2)?N?t+?uN?x+?vN?y=?gh?Z?x?gn2uu2+v2√h1/3(3)?Μ?t+?uΜ?x+?vΜ?y=-gh?Ζ?x-gn2uu2+v2h1/3(2)?Ν?t+?uΝ?x+?vΝ?y=-gh?Ζ?x-gn2uu2+v2h1/3(3)式中:Z為水位,m;h為水深,m;u,v為垂線(xiàn)平均流速在x,y方向的分量,m/s;M,N為單寬流量在x,y方向的分量,M=hu,N=hv,m2/s;n為曼寧糙率系數(shù);g為重力加速度,m/s2.本模型采用有限體積法計(jì)算,離散的詳細(xì)過(guò)程見(jiàn)文獻(xiàn).有限體積法兼有有限差分法物理概念清晰和有限元法適應(yīng)不規(guī)則網(wǎng)格、復(fù)雜邊界情況及計(jì)算精度高等方面的特點(diǎn).有限體積法又稱(chēng)有限控制容積積分法,即把計(jì)算區(qū)域按劃分的網(wǎng)格離散為若干點(diǎn),以這些點(diǎn)為中心,把整個(gè)計(jì)算區(qū)域劃分為若干互相連接但不重疊的控制體.計(jì)算物理量采用交錯(cuò)網(wǎng)格,其優(yōu)越性在于能使水量保持守恒.1.2角形網(wǎng)格單元采用斜對(duì)角笛卡爾法對(duì)實(shí)際水邊界處理后,計(jì)算邊界網(wǎng)格單元較傳統(tǒng)笛卡爾鋸齒法多4種類(lèi)型,共8種基本類(lèi)型(如圖1所示).由圖1可知,邊界網(wǎng)格單元不外乎四邊形和三角形兩種類(lèi)型,計(jì)算時(shí)可通過(guò)網(wǎng)格標(biāo)示數(shù)組判別.基本方程在網(wǎng)格上離散后,邊界三角形網(wǎng)格單元對(duì)應(yīng)的水流參數(shù)(單元形心)與原四邊形網(wǎng)格單元不同,其分布在三角形形心處,而邊界處流速點(diǎn)存儲(chǔ)位置不變.為了求解計(jì)算方法的統(tǒng)一,需對(duì)三角形邊界網(wǎng)格單元進(jìn)行處理,若把三角形邊界網(wǎng)格單元斜邊中點(diǎn)視為一角點(diǎn),則仍可視為四邊形網(wǎng)格,與內(nèi)部計(jì)算網(wǎng)格無(wú)異,可用統(tǒng)一方程求解.由于三角形網(wǎng)格單元水深存儲(chǔ)位置變化,而流速點(diǎn)存儲(chǔ)位置不變,故對(duì)方程組離散參數(shù)影響主要體現(xiàn)在hi,j+1/2,hi+1,j+1/2,hi+1/2,j,hi+1/2,j+1.以hi,j+1/2為例,其值受不同邊界類(lèi)型組合影響,可通過(guò)加權(quán)平均的方法確定各類(lèi)邊界的hi,j+1/2值,基本類(lèi)型有8種,如圖2所示.(a),(b),(d),(e)邊界類(lèi)型有:hi,j+1/2=12(hi?1/2,j+1/2+hi+1/2,j+1/2)hi,j+1/2=12(hi-1/2,j+1/2+hi+1/2,j+1/2)(c),(f)邊界類(lèi)型有:hi,j+1/2=(35hi?1/2,j+1/2+25hi+1/2,j+1/2)hi,j+1/2=(35hi-1/2,j+1/2+25hi+1/2,j+1/2)(c),(f)邊界類(lèi)型有:hi,j+1/2=(25hi?1/2,j+1/2+35hi+1/2,j+1/2)hi,j+1/2=(25hi-1/2,j+1/2+35hi+1/2,j+1/2)三角形網(wǎng)格邊界處其他點(diǎn)處水深如hi+1,j+1/2,hi+1/2,j,hi+1/2,j+1處理與之相同.1.3大陸地區(qū)邊界為檢驗(yàn)采用斜對(duì)角笛卡爾方法改進(jìn)后的數(shù)學(xué)模型較傳統(tǒng)笛卡爾方法優(yōu)越性,本文選取一概化的具有復(fù)雜邊界的分洪區(qū)為計(jì)算對(duì)象,如圖3所示.計(jì)算時(shí),水流邊界條件為進(jìn)口10m寬扒口,分洪流量恒為50m3/s,外邊界為四周封閉圍堤.陸地邊界的處理方法為閉邊界的法向流速為0;初始分洪區(qū)內(nèi)為平底干河床.圖4給出了兩種方法計(jì)算的分洪1000s后區(qū)內(nèi)流場(chǎng)及邊界網(wǎng)格概化示意圖.由圖可知,邊界采用斜對(duì)角笛卡爾方法處理時(shí),其計(jì)算水邊線(xiàn)優(yōu)于鋸齒法的水邊線(xiàn),擬合精度更高;兩種方法計(jì)算得出的流場(chǎng)差別較大,主流線(xiàn)位置發(fā)生改變,其中斜對(duì)角笛卡爾方法流場(chǎng)變化較為平緩,主流線(xiàn)位置距扒口處稍遠(yuǎn).從定性上看,圖3(a)所示流場(chǎng)較圖3(b)的流場(chǎng)更符合實(shí)際情況.結(jié)果表明,斜對(duì)角笛卡爾方法能更好地模擬復(fù)雜邊界水流運(yùn)動(dòng),使邊界流場(chǎng)更加合理.2計(jì)算示例2.1外矛盾分洪區(qū)龍骨湖、沉底湖蓄洪區(qū)位于漢江支流漢北河以北、皂市河以東的湖北省天門(mén)市境內(nèi).蓄滯洪區(qū)總面積51.81萬(wàn)m2,蓄洪水位29.7m(吳淞凍結(jié)高程),蓄洪容積11142萬(wàn)m3.外垸圍堤由漢北河堤(北堤)和皂市河堤(東堤)組成,堤長(zhǎng)度分別為7km和15km,內(nèi)垸沒(méi)有攔洪堤工程設(shè)施,分洪時(shí)洪水按地面等高線(xiàn)漫淹.該分蓄洪區(qū)屬于漢北河防御特大洪水的臨時(shí)應(yīng)急措施,至今未進(jìn)行過(guò)安全轉(zhuǎn)移工程建設(shè),分洪區(qū)示意圖如圖4所示.2.2邊界的處理方法(1)邊界問(wèn)題.邊界條件為:水流邊界條件為扒口分洪時(shí)的外江水位,即分洪水位29.7m,外邊界為分蓄洪民垸的四周?chē)?陸地邊界的處理方法是令閉邊界的法向流速為0.(2)干河床處理.隨著水位的變化,陸地邊界位置不斷變化,對(duì)于水流還沒(méi)有到達(dá)的干河床區(qū)域,需進(jìn)行特殊處理.在分洪區(qū)計(jì)算過(guò)程中,一般作近視簡(jiǎn)化處理.本文采用“凍結(jié)法”處理,根據(jù)計(jì)算水位和湖底高程判斷網(wǎng)格單元是否露出水面,即定義臨界水深Δh=0.005～0.01m,當(dāng)水深h>Δh時(shí),糙率取正常值,反之糙率取一大值(1010量級(jí)).2.3模塊參數(shù)的測(cè)定采用本模型對(duì)龍骨湖、沉底湖蓄洪區(qū)分洪后洪水演進(jìn)進(jìn)行計(jì)算時(shí),由于該分蓄洪區(qū)內(nèi)沒(méi)有實(shí)測(cè)分洪資料,計(jì)算糙率參照相關(guān)文獻(xiàn);以下分析了分蓄洪民垸進(jìn)洪流量、水位、分洪歷時(shí)及流場(chǎng)等變化情況,各特征點(diǎn)位置如圖5所示.(1)分洪動(dòng)態(tài)過(guò)程分析分洪時(shí)分洪口門(mén)處流量過(guò)程見(jiàn)圖6.由圖可知,流量過(guò)程呈現(xiàn)從小變大再到小的過(guò)程,這與實(shí)際分洪規(guī)律相一致,其中最大分洪流量3400m3/s,與設(shè)計(jì)分洪流量相吻合.為檢驗(yàn)?zāi)Ｐ陀?jì)算的水量守恒性,根據(jù)不同瞬時(shí)分洪區(qū)內(nèi)各網(wǎng)格上的水深分布計(jì)算結(jié)果統(tǒng)計(jì)了分洪區(qū)內(nèi)的蓄洪過(guò)程,如圖7所示.同時(shí)根據(jù)圖6的分洪區(qū)進(jìn)洪流量過(guò)程線(xiàn)統(tǒng)計(jì)了分洪區(qū)內(nèi)的進(jìn)洪量過(guò)程.兩者的統(tǒng)計(jì)結(jié)果比較如下:當(dāng)分洪5h時(shí),分洪口門(mén)處進(jìn)洪量為0.45億m3,民院內(nèi)蓄洪量也為0.45億m3;分洪10h,分洪口門(mén)處進(jìn)洪量和民院內(nèi)蓄洪量均為0.895億m3.由此可知,模型在計(jì)算過(guò)程中滿(mǎn)足水量守恒條件.(2)分洪時(shí)段和分洪時(shí)段的特征點(diǎn)蓄洪水位隨分洪時(shí)間的變化圖8給出了A,B,C三特征點(diǎn)蓄洪水位變化過(guò)程,計(jì)算結(jié)果表明:分洪后洪水到達(dá)各特征點(diǎn)的時(shí)間不同,其中距分洪口門(mén)較近的A點(diǎn),洪水在分洪后0.64h到達(dá)該位置,而距口門(mén)較遠(yuǎn)處的B點(diǎn),分洪后6.9h洪水方到達(dá)此位置.各特征點(diǎn)蓄洪水位隨分洪時(shí)間增加緩慢抬升,當(dāng)蓄洪民垸蓄滿(mǎn)后,水位趨于相對(duì)平穩(wěn),分洪停止,分洪歷時(shí)14h.(3)口門(mén)處流場(chǎng)分析分洪時(shí),由于口門(mén)流速通常較大,通常會(huì)給附近建筑帶來(lái)較大影響甚至構(gòu)成威脅,因此必須對(duì)其流場(chǎng)情況作出分析.本文在計(jì)算結(jié)果中選取了蓄滯洪區(qū)在分洪2.8h和11.1h分洪口門(mén)處局部流場(chǎng)進(jìn)行了分析.圖9為分洪口門(mén)分洪2.8h、11.1h瞬時(shí)流場(chǎng),由圖可知,兩時(shí)刻口門(mén)處流場(chǎng)有較大差別,分洪2.8h時(shí),口門(mén)處流量較大,流速也較大,且流場(chǎng)紊亂;分洪11.1h時(shí)分洪流量較小,水位變化緩慢,并趨于穩(wěn)定,流速較小.通過(guò)對(duì)分蓄洪民垸進(jìn)洪流量、水位、分洪歷時(shí)及流場(chǎng)等變化情況分析可知,本文所建的模型能很好地模擬龍骨湖、沉底湖分蓄洪區(qū)分洪后洪水演進(jìn),并為進(jìn)一步的分洪影響評(píng)價(jià)打下基礎(chǔ).3計(jì)算網(wǎng)格、擬合精度差針對(duì)二維水流計(jì)算中復(fù)雜邊界問(wèn)題,本文引進(jìn)了斜對(duì)角笛卡爾方法,該方