淺水方程大時間步長離散方法研究進展

許仁義，趙 磊，王遠見，徐 波

（1.揚州大學 水利科學與工程學院，江蘇 揚州 225009；2.黃河水利科學研究院，河南 鄭州 450003）

在實際的工程項目中，對數(shù)值模擬計算效率的需求是永無止境的。比如上游發(fā)生潰壩，洪峰何時到達下游城市，迅速劃出人口撤離范圍；或者上游發(fā)生化工污染事件，污染物何時會到達下游取水口。這需要流域管理單位根據(jù)模擬軟件的模擬結果迅速做出反應。這時，模擬速度就顯得尤為重要。模擬計算速度的提升方法有兩個方面，一方面是提升硬件水平，在科技日新月異的今天，計算機的主頻提升速度相對放緩，硬件方面可以通過多核來提升計算效率；另一方面是追求好的高效算法，大時間步長格式（Large Time Step，LTS）［1］就是其中一種，它的主要思想是在一個時間步長內完成傳統(tǒng)格式在多個時間步長內完成的工作，這種算法被證明能大幅提升計算效率。本文第一部分從線性標量雙曲線方程求解開始，介紹了大時間步長格式的基本思想和方法，然后將其擴展到非線性雙曲方程組求解，以及二維非結構網(wǎng)格上的大時間步長格式更新方法；第二部分提出了目前大時間步長格式存在的問題以及解決方法，最后對其未來進行了展望。

1 離散方法

1.1 大時間步長格式的引出——標量方程的大時間步長格式

針對線性波動方程（雙曲型偏微分方程）：

式中：u為振幅；a為波的傳播速率，大約為330 m／s。

如果采用如下一階迎風型顯格式離散方法：

那么，需要滿足CFL定理，柯朗數(shù)小于等于1：

如圖1所示，在i單元與i－1單元的界面xi－1／2處產(chǎn)生的波經(jīng)過一個時間步長Δt后到達ε點。采用式（2）離散時，這個ε點需要落在i單元內。

圖1 傳統(tǒng)格式單個時間步長內單波傳輸只限于本單元內

波在時間步長Δt內所行走的距離必然不能超過一個單元的長度，也就是空間步長Δx。

而大時間步長格式是指考慮波在一個時間步長Δt內穿過多個單元的情況，如圖2所示。

圖2 LTS格式單個時間步長內單波傳輸穿過多個單元

此時式（2）的更新方法不再適用。把式（2）改變一個寫法：

式中：ζ為大時間步長格式的更新系數(shù)。

從式（5）可以看出，從n時刻到n＋1時刻，在這個i－1／2處發(fā)出的單波的影響下，變化量是），這個變化量由兩部分組成，其中一個是：

Δu為波強，在傳播過程中是不會變化的。波速a乘以時間Δt得到單波所行進的路程，除以空間步長，結果就是行進路程占整個空間步長的比重，在使用傳統(tǒng)格式時，這個就是柯朗數(shù)，需要滿足小于1。當擴展到大時間步長格式時（見圖2），在界面xi－1／2處產(chǎn)生的單波經(jīng)過一個時間步長Δt后，完整地穿過了i、i＋1單元，最終到達的ε點落在i＋2單元內。由前面的定義，當單波完整地穿過某個單元時，波在單元內所行進的路程與單元的空間步長是相等的。

floor是舍去法求整函數(shù)。比如Cr是2.6，那么floor（Cr）就是2，代入上式得到ζi＋2＝0.6。因此下一時刻的更新方法就變成：

這就是由Leveque［1］提出的大時間步長格式。這種方法的優(yōu)點是：把原來幾個時間步長才能實現(xiàn)的計算工作，在一個時間步長內實現(xiàn)，極大地提升了計算效率。

1.2 非線性雙曲守恒律的大時間步長格式

對于非線性的雙曲守恒律，以淺水方程為例：

式中：h為水深，m；u為流速，m／s；g為重力加速度，取9.8 m／s2。

由標量方程實現(xiàn)大時間步長格式的計算，需要知道波速來求得單波在一個時間步長內穿過了多少單元，所經(jīng)過的單元都要進行更新；需要知道波強以求得這些單元需要怎樣的更新。而這兩個量的求得也有兩種方法，第一種方法是黎曼近似解，就是通常所說的Roe方法［2］。對于相鄰兩個單元先做Roe平均：

由平均值求得雅可比矩陣：

再對雅可比矩陣進行特征化：

求得兩個單元間的界面所發(fā)出波的速度：

求波強：

有了波強和波速，就可以按照式（9）對每個單元進行更新。

第二種方法是黎曼精確解，這種方法需要迭代。對于相鄰的兩個單元，有著不同的初始變量，可以理解為一個局部間斷，那么這就是一個黎曼問題。非線性雙曲守恒系統(tǒng)的黎曼問題的解有4種形式［3］，它們是稀疏波與激波的4種排列組合：①左右兩個都是激波；②左稀疏波右激波；③左激波右稀疏波；④兩側都為稀疏波。問題求解的關鍵是求得中間狀態(tài)。如果是激波，則滿足條件：

式中：S為激波的速度。F和U見（10）式，下標S表示中間狀態(tài)，下標L表示左邊，下標R表示右邊。

如果左邊是稀疏波，則：

如果右邊是稀疏波，則：

通過式（17）～式（19）可以求出中間狀態(tài)，與兩側已知的左右狀態(tài)的差值即為波強，而波速也可以通過上面的關系式求出。

1.3 含底坡源項的大時間步長格式

含底坡源項的淺水方程，若采用簡單的中心格式則會產(chǎn)生虛假流動，傳統(tǒng)格式為解決這個問題采用的是和諧（well?balanced）格式［4－6］。大時間步長格式對此也有兩種方法，第一種方法是黎曼近似解，對底坡源項進行特征分解［7］，把方程（10）變成：

式中：z為河底高程。

此時雅可比矩陣的右矩陣為

此時雖然由前面的兩個單波變成3個，但是其中一個波速為0，另外兩個波速不變：

波強為

那么：

第二種方法同樣是黎曼精確解。采用Godonov方法離散時，兩個相鄰單元的物理狀態(tài)是不連續(xù)的，同樣在非齊次的淺水方程中，底坡也是不連續(xù)的，像樓梯的臺階一樣，我們稱之為階梯黎曼問題（Step Riemann Problem，SRP）。相比于齊次黎曼問題，階梯黎曼問題復雜許多，因為在階梯位置出現(xiàn)了一個靜態(tài)的激波，所以待求的中間狀態(tài)變成了兩個，遵循的關系依然是式（18）～式（20）。

1.4 二維大時間步長格式

將一維格式拓展到結構網(wǎng)格的二維大時間步長格式，采用Strane維分裂的方法［8］：當時間步長n是奇數(shù)時，x方向上前進半個時間步長，然后y方向上前進一個時間步長，最后回到x方向上再前進半個時間步長；當時間步長n是偶數(shù)時，y方向上前進半個時間步長，然后x方向上前進一個時間步長，最后回到y(tǒng)方向上再前進半個時間步長。

如果是二維非結構網(wǎng)格，問題就會變得復雜［9－10］，還是按照前面的方法，先求出波速和波強。對于二維淺水方程：

式中：v為橫向流速；＝icosθ＋jsinθ為單元邊的外法線方向。

兩個相鄰的單元做Roe平均得到：

于是得到雅可比矩陣：

對其進行特征化：

左矩陣為

右矩陣為

波速為

波強為

于是得到：

得到波速和波強的表達式后，下面來看非結構網(wǎng)格上的大時間步長格式是如何實現(xiàn)的。

如圖3所示的非結構三角形網(wǎng)格，在三角形單元Cw上一邊w產(chǎn)生的單波，如果是傳統(tǒng)格式，只會影響Cw單元本身，而在大時間步長格式里還要擴展到Cw單元的相鄰單元C1和C2。與一維大時間步長格式不同的是，在波的傳播方向上，一維格式的相鄰單元只有一個，而二維三角形網(wǎng)格有兩個，這就需要對波強進行分配。

圖3 二維非結構網(wǎng)格上的LTS

2 存在的問題

2.1 非線性問題

大時間步長格式從線性標量方程擴展到非線性雙曲守恒律遇到的第一個問題是非線性的問題，在大時間步長格式中，每個單元需要經(jīng)歷不同界面的波，如圖4所示，如果波是線性的或者非線性比較弱，那么只需要將這些波的波強簡單累加就可以；如果非線性比較強，那么不同波相遇后產(chǎn)生新波的波速會發(fā)生變化，如圖5所示。空氣動力學的歐拉方程中，是需要考慮這種非線性疊加的［11－12］，而在淺水方程中，考慮這種非線性帶來的變化幾乎可以忽略［13］。

圖4 多個單波同時穿過一個單元

圖5 非線性單波疊加

2.2 稀疏波的分解

采用Godonov型格式進行離散時，每兩個單元間的界面上是一個局部間斷，就是一個黎曼問題。對非線性雙曲系統(tǒng)，這個黎曼問題會產(chǎn)生稀疏波和激波。我們通常使用的Roe［2］、Osher［14］、HLL［15］等格式是指對這個黎曼問題的近似解。在這些傳統(tǒng)格式中，忽略了稀疏波的存在，將其視為激波。因為時間步長不夠大，所以稀疏波的寬度擴張不是很嚴重，如圖6中的紅色部分。但是在大時間步長格式中，這個問題就凸顯出來。如果仍然采用傳統(tǒng)方法不做任何處理，就會出現(xiàn)稀疏波分裂的情況（如圖7所示）。

圖6 時間步長增大稀疏波波頭波尾距離變大

圖7 稀疏波斷裂

Morales－Hernández等［16］在黎曼近似解中強制對單波進行分解，也取得了不錯的效果。另一種方法就是采用黎曼精確解。這種方法在傳統(tǒng)的格式中沒有被廣泛使用的原因是在每個局部間斷需要做Newton迭代來取得精確解，而且迭代的步數(shù)無法控制，可能造成迭代步數(shù)過多甚至發(fā)散。而這種方法最大的好處就是可以完美地求解出稀疏波，得到波頭和波尾的速度，使得稀疏波的分解變得簡單。實踐證明，大時間步長格式中迭代步數(shù)一般為3步左右，而對于柯朗數(shù)動輒幾十的大時間步長格式來說，這個計算成本完全可以抵消［17］。

2.3 底坡源項

淺水方程有底坡和糙率兩個源項，使得原來的齊次雙曲偏微分方程組變成了非齊次，而大時間步長格式解決這個問題有兩大類方法。第一大類是黎曼近似解，因為和諧格式是用物理通量平衡源項，所以這種方法是不能采用傳統(tǒng)和諧格式的，而此方法中，物理通量沒有被計算，那只有對底坡源項進行特征分解，將其融合到波強中，隨著單波一起穿過多個單元［18］。

第二大類是黎曼精確解，這類方法比較復雜。首先求解的理論就分為兩種，一種是基于動量和質量守恒的方法［19］，一種是基于能量質量守恒的求解方法［20］；其次是對底坡源項間斷的處理方法，一種認為這是一個垂直間斷［21－22］，一種認為這是一個單調連續(xù)的變量［23－24］（寬度趨向于0）。階梯黎曼問題在很多情況下不能得到唯一解，只有通過附加其他條件才能使得解唯一［25－26］。因為求解的出發(fā)點不同，所以這些方法求得的結果不能統(tǒng)一。

2.4 振 蕩

振蕩問題可以說是大時間步長格式的頑疾，采用隨機選取法（Random Choice Method，RCM）只能在一定程度上抑制振蕩［27］，隨著柯朗數(shù)的增大，這種方法也逐漸失去了效用。在計算機程序中，隨機數(shù)不是真正的隨機數(shù)，而是固定的。Wang等［28］用范德科皮特（Van der Corput）序列來制造隨機數(shù)，效果沒有什么變化。如果采用濾波法來抑制振蕩，這顯然需要預先知道振蕩發(fā)生的位置和幅度，在工程中不太可能實現(xiàn)。所以，到目前為止對這個問題還沒有一個很好的解決方法。

3 總結與展望

大時間步長格式從提出至今已有40多a時間，最近幾年仍然有很多成果，展現(xiàn)了其強大的生命力。從另一方面來看，它仍然有問題沒有得到解決。其基本思想就是通過單個時間步長內單波穿過多個單元，在一個時間步長內完成了傳統(tǒng)格式多個時間步長的計算，通過這種方法可極大地提高計算效率，這在線性波動方程里是顯而易見的，但是到了非線性的方程組里就會出現(xiàn)一系列問題，其中最重要的問題就是振蕩，這種振蕩與傳統(tǒng)格式的振蕩在現(xiàn)象上有所區(qū)別，不是發(fā)生在激波附近，而是分布在整個平臺。其解決方法也不同于傳統(tǒng)格式，這是個頑疾，如果這個問題得到解決，大時間步長格式將會得到更廣闊的應用。