陳善群，廖 斌，李海峰

(1.安徽工程大學(xué)建筑工程學(xué)院，安徽 蕪湖 241000;2.中廣核風(fēng)力發(fā)電有限公司華南分公司，廣東深圳 518031)

1 研究背景

在以海岸、河口、湖泊、大型水庫(kù)等廣闊水域地區(qū)為模型，進(jìn)行流體數(shù)值模擬時(shí)，由于水平尺度遠(yuǎn)大于垂向尺度，水力參數(shù)(如流速、水深等)在垂直方向變化要明顯小于水平方向的變化，其流態(tài)可用沿水深的平均流動(dòng)量來(lái)表示。對(duì)于這些模型的數(shù)值模擬，可采用平面二維水動(dòng)力數(shù)值模擬技術(shù)。二維水動(dòng)力數(shù)值模擬技術(shù)中，使用的控制方程是將三維流動(dòng)基本方程沿水深積分平均得到，又稱為二維淺水方程。

對(duì)于二維淺水方程的求解，國(guó)內(nèi)外學(xué)者已做了大量研究，目前主要的研究集中在2個(gè)方面:基于非結(jié)構(gòu)網(wǎng)格的空間離散［1－2］和高性能的計(jì)算格式［3－4］。近年來(lái)興起的無(wú)網(wǎng)格方法［5］具有無(wú)需劃分網(wǎng)格、克服場(chǎng)函數(shù)局部化近似所引起的誤差、僅需對(duì)邊界條件進(jìn)行描述、無(wú)需尋求光滑梯度場(chǎng)的后處理、適用于涉及大變形和需要?jiǎng)討B(tài)調(diào)整的各類應(yīng)用問(wèn)題、適合進(jìn)行自適應(yīng)分析，以及求解精度較高等優(yōu)點(diǎn)，已在固體力學(xué)數(shù)值計(jì)算領(lǐng)域得到成功的運(yùn)用。在計(jì)算水力學(xué)，特別在求解二維淺水方程方面，運(yùn)用無(wú)網(wǎng)格方法在國(guó)內(nèi)還幾乎是一片空白。

本文擬將最小二乘無(wú)網(wǎng)格有限差分方法(簡(jiǎn)稱MLSFD方法)［6］運(yùn)用于求解二維淺水方程，通過(guò)算例驗(yàn)證無(wú)網(wǎng)格方法在求解二維淺水方程中的可行性。

2 無(wú)網(wǎng)格有限差分方法的相關(guān)概念

2.1 支撐域［5］

以二維空間為例，如圖1所示，將求解域Ω用n個(gè)節(jié)點(diǎn)離散，域中所有圓點(diǎn)即為相互獨(dú)立的計(jì)算節(jié)點(diǎn);在求解域計(jì)算節(jié)點(diǎn)中任找一個(gè)中心點(diǎn)(同心圓點(diǎn))，則此中心點(diǎn)的函數(shù)值可由它周圍臨近區(qū)域的計(jì)算節(jié)點(diǎn)通過(guò)一定的代數(shù)組合來(lái)逼近，我們稱這種臨近區(qū)域Ω為中心點(diǎn)的支撐域。

圖1 節(jié)點(diǎn)支撐域Fig.1 Support domain of nodes

無(wú)網(wǎng)格方法的支撐域理論上可以是各種平面圖形，但為了便于計(jì)算，在具體選擇時(shí)一般取大小合適的圓形或矩形區(qū)域。圖1中有2個(gè)支撐域，以圓形支撐域?yàn)槔膱A點(diǎn)為支撐域中心點(diǎn)，實(shí)心黑點(diǎn)為支撐域內(nèi)點(diǎn)，這兩者一同組成“點(diǎn)云”參與中心點(diǎn)函數(shù)值的計(jì)算;而支撐域外點(diǎn)(空心圓點(diǎn))，我們認(rèn)為其不參與中心點(diǎn)函數(shù)值的計(jì)算，對(duì)中心點(diǎn)的影響為0。

2.2 支撐域內(nèi)點(diǎn)的選取原則

支撐域范圍確定后，不是所有支撐域內(nèi)的點(diǎn)都參與中心點(diǎn)的擬合計(jì)算。本文采用2種方法對(duì)支撐域內(nèi)的點(diǎn)進(jìn)行了選?。?］:

(1)擬定參與計(jì)算的點(diǎn)為n個(gè)，取支撐域內(nèi)距離中心點(diǎn)最近的n個(gè)點(diǎn)即可。這種方法相當(dāng)簡(jiǎn)單，耗機(jī)時(shí)較少;

(2)在前一種方法的基礎(chǔ)上，對(duì)這n個(gè)點(diǎn)再進(jìn)行一輪篩選。所選擇點(diǎn)需盡量均勻分布在中心點(diǎn)周圍(圖2(a))，應(yīng)避免在某個(gè)大扇形區(qū)域內(nèi)無(wú)點(diǎn)(圖2(b))，也要避免過(guò)小角度區(qū)域內(nèi)有多個(gè)點(diǎn)(圖2(c))，取扇形角度 在40°～120°之間。保證參與中心點(diǎn)計(jì)算的節(jié)點(diǎn)數(shù)為4～9個(gè)。在2個(gè)節(jié)點(diǎn)與中心點(diǎn)連接后夾角過(guò)小，需要舍棄1個(gè)時(shí)，保留離中心點(diǎn)近的節(jié)點(diǎn)(圖2(c))，保留節(jié)點(diǎn)1，去掉節(jié)點(diǎn)2。這種方法相對(duì)前一種方法要耗機(jī)時(shí)，但計(jì)算點(diǎn)分布均勻。計(jì)算節(jié)點(diǎn)選取數(shù)量因中心點(diǎn)支撐域中離散點(diǎn)分布疏密的不同，為一不確定數(shù)。

圖2 支撐域內(nèi)點(diǎn)選取示意圖Fig.2 Selection of interior nodes in support domain

2.3 最小二乘法的運(yùn)用

無(wú)網(wǎng)格有限差分法計(jì)算導(dǎo)數(shù)近似值時(shí)會(huì)遇到這種情況:由于參與計(jì)算的支撐域內(nèi)點(diǎn)相互之間非?？拷斐傻南禂?shù)矩陣奇異和病態(tài)。這種問(wèn)題并不容易解決，因?yàn)槲覀儾⒉涣私庵斡騼?nèi)的節(jié)點(diǎn)相互間的空間分布對(duì)系數(shù)矩陣造成的影響，除非這些節(jié)點(diǎn)在一條直線上。為了讓數(shù)值計(jì)算進(jìn)行下去，就得檢驗(yàn)和確保系數(shù)矩陣在計(jì)算區(qū)域內(nèi)的每個(gè)節(jié)點(diǎn)處都是良態(tài)的、可逆的，然而這種做法會(huì)大大增加計(jì)算機(jī)時(shí)。

為了避開(kāi)繁瑣的矩陣奇異、病態(tài)的檢驗(yàn)過(guò)程，此處選用了最小二乘技術(shù)對(duì)矢量做最優(yōu)化近似，具體實(shí)現(xiàn)就是選取足夠多的支撐域內(nèi)節(jié)點(diǎn)(大于要求解的導(dǎo)數(shù)數(shù)量)，組成方程組進(jìn)行計(jì)算，從而避開(kāi)系數(shù)矩陣的奇異問(wèn)題。

3 控制方程

本文選取控制方程為二維淺水方程非守恒形式，忽略柯氏力和風(fēng)對(duì)模型的影響，最終形式為

式中:u，v分別為流速在x，y方向的分量;h0為靜水時(shí)的水深;ξ為自由水面在豎直方向的位移;τbx，τby分別為床面阻力在x，y方向的分量。

4 控制方程的離散

MLSFD方法的離散形式［8］為:

將式(2)代入式(1)的空間導(dǎo)數(shù)中，時(shí)間上仍采用向前差分格式，得到二維淺水方程的MLSFD離散式:

5 圓形潰壩模型

數(shù)值計(jì)算中，常用圓形潰壩模型來(lái)檢驗(yàn)算法在處理對(duì)稱間斷流方面的能力。本文將使用MLSFD(Meshfree Least-Squares-based Finite Difference Method)方法對(duì)該模型進(jìn)行計(jì)算。

5.1 計(jì)算模型具體設(shè)置及初始條件

模型具體設(shè)置及初始條件:計(jì)算區(qū)域?yàn)橐婚L(zhǎng)和寬都為50 m的正方形區(qū)域，區(qū)域中心設(shè)置一半徑11 m，水深10 m的圓形水壩，其他區(qū)域水位為1 m，在某一時(shí)刻圓形水壩突然潰決。圖3為圓形潰壩模型初始時(shí)刻水位示意圖。為提高計(jì)算精度，時(shí)間步長(zhǎng)取0.001 s。

圖3 圓形潰壩初始水位Fig.3 Initial water level for the circular dam-break model

5.2 計(jì)算模型網(wǎng)格劃分及邊界條件

如圖4所示，計(jì)算區(qū)域中心采用無(wú)網(wǎng)格隨機(jī)布點(diǎn)，計(jì)算節(jié)點(diǎn)總數(shù)為11 823個(gè)。計(jì)算模型四周邊界條件為:x 方向?u/?x=0，?v/?x=0，?ξ/?x=0;y 方向?u/?y=0，?v/?y=0，?ξ/?y=0。

圖4 圓形潰壩模型網(wǎng)格劃分Fig.4 Mesh generation for the circular dam-break model

5.3 計(jì)算模型邊界處理

為了更好地處理邊界條件，我們?cè)跓o(wú)網(wǎng)格的四周，向內(nèi)分別布置了3層直角網(wǎng)格。需要注意的是，這種特殊的網(wǎng)格布置，僅僅是用來(lái)處理邊界條件的。除最外層邊界上的點(diǎn)，其他所有點(diǎn)包括直角網(wǎng)格上的節(jié)點(diǎn)，仍采用MLSFD格式計(jì)算。如圖5所示，在邊界上的w點(diǎn)的函數(shù)值，根據(jù)邊界條件，可以由內(nèi)層的w+1，w+2點(diǎn)的函數(shù)值來(lái)求得。令?表示節(jié)點(diǎn)處的函數(shù)值，則邊界條件為??/?xi=0，通過(guò)一維泰勒級(jí)數(shù)在這3點(diǎn)上的展開(kāi)式，可以求得w點(diǎn)的有限差分格式:

圖5 圓形潰壩模型邊界處理Fig.5 Boundary treatment for the circular dam-break model

5.4 計(jì)算結(jié)果分析

計(jì)算時(shí)選取支撐域內(nèi)離中心點(diǎn)距離最近的17個(gè)點(diǎn)，將初始值和邊界條件處理格式代入式(3)，在每一時(shí)間步迭代求解，得出0.69 s時(shí)圓形潰壩模型計(jì)算模擬結(jié)果(圖6、圖7、圖8)。圖6中所示的圓形潰壩0.69 s時(shí)的水位呈尖錐形，頂部水位最高。由于是突然潰決，上面潰決的水體對(duì)下面的水體有向下的沖擊作用，這樣就在水位尖錐底部周圍形成了一個(gè)環(huán)形的淺水坑。圖7中所示的圓形潰壩0.69 s時(shí)的水位等值線呈圓環(huán)狀，分布比較均勻，其中圓環(huán)中心處的水位接近10.0 m，邊緣處水位接近1.0 m。圖6、圖7所示的數(shù)值模擬結(jié)果與 Mingham［9］采用的極坐標(biāo)網(wǎng)格計(jì)算的結(jié)果相當(dāng)吻合。圖8中所示的速度矢量均勻發(fā)散，中心處速度最小，邊緣處速度最大，說(shuō)明圓形水體是從邊緣往中心處潰決，這與實(shí)際情況相符。

圖6 圓形潰壩0.69 s時(shí)的水位Fig.6 Water level for the circular dam break at t=0.69 s

圖7 圓形潰壩0.69 s時(shí)的水位等值線Fig.7 Contour plot of water level for the dam break at t=0.69 s

圖8 圓形潰壩0.69 s時(shí)的速度矢量Fig.8 Velocity vectors for the dam break at t=0.69 s

6 結(jié)論

本文將最小二乘無(wú)網(wǎng)格有限差分(MLSFD)方法運(yùn)用于求解二維淺水方程，利用最小二乘無(wú)網(wǎng)格有限差分離散式對(duì)二維方程進(jìn)行了離散求解。在計(jì)算圓形潰壩這一典型數(shù)值算例時(shí)，對(duì)計(jì)算區(qū)域進(jìn)行離散布點(diǎn)，利用泰勒離散式處理邊界條件，得出了0.69 s時(shí)圓形潰壩的水位圖、水位等值線圖以及速度矢量圖。并對(duì)數(shù)值模擬結(jié)果進(jìn)行分析比較，驗(yàn)證了該方法在計(jì)算二維淺水流動(dòng)的數(shù)值模擬方面有著較高的精確度，進(jìn)而說(shuō)明無(wú)網(wǎng)格方法運(yùn)用于求解淺水方程是可行的。

［1］楊 彬，汪德?tīng)?非結(jié)構(gòu)網(wǎng)格上淺水方程的LU-SGS隱式算法［J］.河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2008，36(4):483 －487.(YANG Bin，WANG De-guan.LU-SGS Scheme for Shallow Water Equation Computation on Unstructured Grids［J］.Journal of Hohai University(Natural Sciences)，2008，36(4):483 －487.(in Chinese))

［2］盧康明，李光熾.非結(jié)構(gòu)網(wǎng)格淺水方程隱式解法［J］.水動(dòng)力學(xué)研究與進(jìn)展(A輯)，2010，25(2):247－253.(LU Kang-ming，LI Guang-chi.An Implicit Method for Shallow Water Equations on Unstructured Grids［J］.Chinese Journal of Hydrodynamics，2010，25(2):247 － 253.(in Chinese))

［3］潘存鴻.三角形網(wǎng)格下求解二維淺水方程的和諧Godunov格式［J］.水科學(xué)進(jìn)展，2007，18(2):204 －209.(PAN Cun-hong.Well-balanced Godunov-type Scheme for 2D Shallow Water Flow with Triangle Mesh［J］.Advances in Water Science，2007，18(2):204 － 209.(in Chinese))

［4］潘存鴻，徐 昆.三角形網(wǎng)格下求解二維淺水方程的KFVS格式［J］.水利學(xué)報(bào)，2006，37(7):858－864.(PAN Cun-hong，XU Kun.Kinetic Flux Vector Splitting Scheme for Solving 2D Shallow Water Equations with Triangular Mesh［J］.Journal of Hydraulic Engineering，2006，37(7):858 －864.(in Chinese))

［5］張 雄，劉 巖.無(wú)網(wǎng)格法［M］.北京:清華大學(xué)出版社，2004.(ZHANG Xiong，LIU Yan.Meshless Methods［M］.Beijing:Tsinghua University Press，2004.(in Chinese))

［6］DING H，SHU C.Simulation of Incompressible Viscous Flows Past a Circular Cylinder by Hybrid FD Scheme and Meshless Least Square-based Finite Difference Method［J］.Computer Methods in Applied Mechanics and Engineering，2004，193:727 －744.

［7］江興賢，陳紅全，Euler方程無(wú)網(wǎng)格算法及布點(diǎn)技術(shù)［J］.南京航空航天大學(xué)學(xué)報(bào)，2004，36(2):174 －178.(JIANG Xing-xian，CHEN Hong-quan.Gridless Method for Euler Equations and Distributing Point Technique［J］.Journal of Nanjing University of Aeronautics ＆Astronautics，2004，36(2):174 －178.(in Chinese))

［8］DING H，SHU C.Development of Least-Square-Based Two-Dimensional Finite-Difference Schemes and Their Application to Simulate Natural Convection in a Cavity［J］.Computers＆ Fluids，2004，(33):137 －154.

［9］MINGHAM C G，CAUSON D M.High-Resolution Finite-Volume Method for Shallow Water Flows［J］.Journal of Hydraulic Engineering，1998，124(6):605－614.