用最小二乘無網格有限差分方法求解二維淺水方程

陳善群，廖 斌，李海峰

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

1 研究背景

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

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

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

2 無網格有限差分方法的相關概念

2.1 支撐域［5］

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

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

無網格方法的支撐域理論上可以是各種平面圖形，但為了便于計算，在具體選擇時一般取大小合適的圓形或矩形區(qū)域。圖1中有2個支撐域，以圓形支撐域為例，同心圓點為支撐域中心點，實心黑點為支撐域內點，這兩者一同組成“點云”參與中心點函數(shù)值的計算;而支撐域外點(空心圓點)，我們認為其不參與中心點函數(shù)值的計算，對中心點的影響為0。

2.2 支撐域內點的選取原則

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

(1)擬定參與計算的點為n個，取支撐域內距離中心點最近的n個點即可。這種方法相當簡單，耗機時較少;

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

圖2 支撐域內點選取示意圖Fig.2 Selection of interior nodes in support domain

2.3 最小二乘法的運用

無網格有限差分法計算導數(shù)近似值時會遇到這種情況:由于參與計算的支撐域內點相互之間非?？拷?，造成的系數(shù)矩陣奇異和病態(tài)。這種問題并不容易解決，因為我們并不了解支撐域內的節(jié)點相互間的空間分布對系數(shù)矩陣造成的影響，除非這些節(jié)點在一條直線上。為了讓數(shù)值計算進行下去，就得檢驗和確保系數(shù)矩陣在計算區(qū)域內的每個節(jié)點處都是良態(tài)的、可逆的，然而這種做法會大大增加計算機時。

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

3 控制方程

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

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

4 控制方程的離散

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

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

5 圓形潰壩模型

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

5.1 計算模型具體設置及初始條件

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

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

5.2 計算模型網格劃分及邊界條件

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

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

5.3 計算模型邊界處理

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

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

5.4 計算結果分析

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

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

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

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

6 結論

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

［1］楊 彬，汪德爟.非結構網格上淺水方程的LU-SGS隱式算法［J］.河海大學學報(自然科學版)，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］盧康明，李光熾.非結構網格淺水方程隱式解法［J］.水動力學研究與進展(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］潘存鴻.三角形網格下求解二維淺水方程的和諧Godunov格式［J］.水科學進展，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］潘存鴻，徐 昆.三角形網格下求解二維淺水方程的KFVS格式［J］.水利學報，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］張 雄，劉 巖.無網格法［M］.北京:清華大學出版社，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方程無網格算法及布點技術［J］.南京航空航天大學學報，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.