許仁義，趙 磊，王遠(yuǎn)見，徐 波

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

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

1 離散方法

1.1 大時間步長格式的引出——標(biāo)量方程的大時間步長格式

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

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

如果采用如下一階迎風(fēng)型顯格式離散方法：

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

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

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

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

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

圖2 LTS格式單個時間步長內(nèi)單波傳輸穿過多個單元

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

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

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

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

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

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

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

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

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

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

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

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

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

求波強：

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

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

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

如果左邊是稀疏波，則：

如果右邊是稀疏波，則：

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

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

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

式中：z為河底高程。

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

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

波強為

那么：

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

1.4 二維大時間步長格式

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

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

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

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

于是得到雅可比矩陣：

對其進行特征化：

左矩陣為

右矩陣為

波速為

波強為

于是得到：

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

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

圖3 二維非結(jié)構(gòu)網(wǎng)格上的LTS

2 存在的問題

2.1 非線性問題

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

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

圖5 非線性單波疊加

2.2 稀疏波的分解

采用Godonov型格式進行離散時，每兩個單元間的界面上是一個局部間斷，就是一個黎曼問題。對非線性雙曲系統(tǒng)，這個黎曼問題會產(chǎn)生稀疏波和激波。我們通常使用的Roe［2］、Osher［14］、HLL［15］等格式是指對這個黎曼問題的近似解。在這些傳統(tǒng)格式中，忽略了稀疏波的存在，將其視為激波。因為時間步長不夠大，所以稀疏波的寬度擴張不是很嚴(yán)重，如圖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］。

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

2.4 振 蕩

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

3 總結(jié)與展望

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