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

陳善群，廖 斌，李海峰

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

1 研究背景

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

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

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

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

2.1 支撐域［5］

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

圖1 節點支撐域Fig.1 Support domain of nodes

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

2.2 支撐域內點的選取原則

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

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

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

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

2.3 最小二乘法的運用

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

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

3 控制方程

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

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

4 控制方程的離散

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

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

5 圓形潰壩模型

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

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

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

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

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

如圖4所示，計算區域中心采用無網格隨機布點，計算節點總數為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層直角網格。需要注意的是，這種特殊的網格布置，僅僅是用來處理邊界條件的。除最外層邊界上的點，其他所有點包括直角網格上的節點，仍采用MLSFD格式計算。如圖5所示，在邊界上的w點的函數值，根據邊界條件，可以由內層的w+1，w+2點的函數值來求得。令?表示節點處的函數值，則邊界條件為??/?xi=0，通過一維泰勒級數在這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時的水位呈尖錐形，頂部水位最高。由于是突然潰決，上面潰決的水體對下面的水體有向下的沖擊作用，這樣就在水位尖錐底部周圍形成了一個環形的淺水坑。圖7中所示的圓形潰壩0.69 s時的水位等值線呈圓環狀，分布比較均勻，其中圓環中心處的水位接近10.0 m，邊緣處水位接近1.0 m。圖6、圖7所示的數值模擬結果與 Mingham［9］采用的極坐標網格計算的結果相當吻合。圖8中所示的速度矢量均勻發散，中心處速度最小，邊緣處速度最大，說明圓形水體是從邊緣往中心處潰決，這與實際情況相符。

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

［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.