波動數值模擬的穩定性

廖振鵬，謝志南

(1.中國地震局工程力學研究所，黑龍江哈爾濱 150080;2.哈爾濱工程大學船舶學院，黑龍江哈爾濱 150001)

波動數值模擬是力學、電磁學、聲學、地球物理和多個工程學科共同關注的領域，研究者在基礎和應用研究方面皆取得了豐碩成果.但是，邊界引入的失穩問題尚未徹底解決，特別是由邊界引入的局部失穩的機理和消除方法，即使在線性范圍內亦待進一步研究.這里“邊界”泛指人工邊界、物理邊界或不同介質的分界面，“局部失穩”則指邊界及其鄰近空間區域在有限時間內由數值解格式不當所引發的誤差放大.這一放大可能導致數值解災難性發散而終止計算，亦可能僅在數值解中引入附加誤差.例如，模擬輻射條件的人工邊界格式與相鄰內域節點的格式匹配不當就可能引發局部失穩［1-2］.迄今為止實用的穩定性分析方法為von Neumann分析［3-4］和譜分析［5］，前者從無限模型(含周期性邊界模型)導出，未考慮非周期邊界影響;后者從有限模型導出，在一定程度上考慮了邊界影響，但未考慮邊界可能引發的局部失穩.科技工作者通常采用經驗和半經驗的方法處理穩定問題［1-2］，數學家則對偏微分方程初邊值問題數值解的穩定性作了深入的研究［6-11］.但是，數學家所采用的簡化模型與各專業領域涉及的真實模型之間存在差距，同時，其分析方法與實際應用所要求的可操作性之間亦有差距.因此，應用后者提供的有益思想研究真實模型中發生的局部失穩問題，并由此發展實用的穩定性分析方法值得注意.本文試圖在線性范圍內考察和梳理穩定性的基本概念和主要研究方法，以利于進一步發展波動數值模擬技術，特別是解決人工邊界引發的局部失穩問題.

1 強穩定性和弱穩定性

從收斂性和Lax等價定理出發，依據收斂方式的不同將Lax穩定性區分為強穩定性和弱穩定性.Lax等價定理是分析線性波動數值模擬問題的合理切入點［12］.該定理可表述為:如果連續模型的初邊值問題適定和數值解格式相容，則收斂性等價于Lax穩定性.Lax穩定性涵義如下:就線性波動數值模擬的理論分析而言，數值解可寫成

式中:矩陣A為數值解算子，取決于內域場方程和邊界條件的時空離散格式，因此，數值解方案的所有性質由A確定.A的階數J與空間步距h相關，J階向量vn和v0分別表示t=nk和t=0時刻的數值解，n為時步數，k為時間步距.Lax穩定性由下式定義:式中:T為計算時間長度，C為與k，h和T無關的常數，k0為某一正數，lnt表示取整.式(2)可表述為A的冪的范數一致有界.因此，穩定性分析是一個代數問題.穩定性分析的目的為確定條件式(2)對k和h取值所施加的限制，即k和h的取值應在k-h平面第一象限的某一區域之上.此區域稱為穩定區，考慮到收斂性，穩定區以k-h平面坐標原點為聚點［11］.式(2)作為一個代數條件與時、空步距趨于零的方式無關，但是，從數值解收斂的角度分析，k和h趨于零的不同方式將導致條件式(2)本質上的差異.由此出發，這里依據k和h趨于零的2種不同方式將Lax穩定性區分為強穩定性和弱穩定性.

1.1 強穩定性

設在網格節點上連續模型的解為u(t)，強穩定性可以通過直接考察數值解vn對連續模型解un=u(nk)的收斂性引入.收斂性由下式定義:

在波動數值解穩定性論著中通常引入參數Δτ=k/h，則穩定區可用參數Δτ和k描述.在數值解格式僅涉及參數Δτ的情況下，由穩定條件式(2)導出的穩定區簡化為允許Δτ取值的實數區間，例如

式中:Δτ*為穩定區的上界.對于穩定區內任一Δτ，h取決于k，故式(2)應在時、空步距同步趨于零時成立.數值解格式的強穩定性定義為在k和h同步趨于零的條件下式(2)成立.因此，數值解格式是否具有強穩定性，只需建立如式(4)所示類型的穩定準則，并檢查Δτ=k/h是否滿足此準則.由于k和h同步趨于零，當k→0時J→∞，即A的階數隨k→0而趨于無窮.由此可知，強穩定性涉及對一簇階數變化的矩陣A的分析.對于模擬真實世界的模型，如何依據強穩定性確定穩定區，就作者所知，在一些重要情況下實用的分析方法尚待完善.

1.2 弱穩定性

弱穩定性可以通過分別考察數值解vn的空間和時間收斂性引入.首先對連續模型作空間離散得離散解w(t)所滿足的常微分方程組.就穩定性分析而言，此方程組可以寫成

式中:J階矩陣B為常微分算子.將式(5)作時間離散可導得式(1).于是，式(1)對于連續模型解的收斂問題簡化為相互獨立的2個問題:

1)式(5)的初值問題的收斂性:

2)時空離散模型的解vn對wn=w(nk)的收斂性.在保證前一收斂性的前提下，數值解格式的弱穩定性定義為:在給定h的條件下當k趨于零時式(2)成立.因為h已給定，A的階數J與k無關，故弱穩定性分析僅涉及階數不變的一簇算子A.依據弱穩定性確定穩定區已有成熟的譜分析方法.

就保證數值解收斂而言，強、弱穩定性要求并無差別.兩者的區別在于:1)前者直接分析時空離散格式，故可用于分析一般數值解格式可能出現的所有失穩現象，特別是局部失穩現象;后者則限用于分析離散模型中時、空變量相互獨立的情形，特別是適用于波動的有限元模擬.2)就分析方法的可操作性而言，后者易于前者，且可利用有限模型的經典譜判據建立穩定準則，而前者在一些情況下建立穩定準則的實用方法尚待研究.3)就數值試驗而言，強穩定性分析提供了穩定區間的界限，便于選取時間步距k;而弱穩定性保證Δτ→0時式(2)成立，為在數值模擬中通過逐步減小時間步距獲得穩定的數值結果提供了依據.雖然弱穩定性難以為k的取值提供依據，這一缺點并不會成為數值試驗的障礙，可以參考von Neumann分析結果或CEL條件試選k的實驗數據［11］.下面進一步闡明這兩類穩定性分析的基本概念及其演進.

2 強穩定性分析

回顧強穩定性譜分析和正則模態分析的主要進展，著重評論 Godunov和 Ryabenkii的開創性工作［6］，特別是揭示在這項工作中隱含的對解決局部失穩問題有價值的思想.

2.1 譜分析的進展及其局限性

有限模型的譜分析建立在A的特征值和特征向量概念之上.A的全部特征值λ1，λ2…，λJ的集合{λi}稱為A的譜，與譜{λi}對應的特征向量的集合記為{vi}.直接用A的譜給出的穩定性判據稱為譜判據.為了獲得保證式(2)成立的譜判據可利用任一方陣的約當變換:

式中:Q為J階非奇異陣，Λ為A的約當標準型.假定

式中:常數C與k和h無關，則得如下譜判據:1)對任意A，穩定的必要條件為

穩定的充分條件為

式中:譜半徑 ρ(A)=max|λi|.

2)當A非缺損(non-defective)時，即當A有J個線性無關的特征向量時，穩定的充要條件為式(9).

3)當A缺損(defective)時，穩定的充要條件為式(9)加上如下條件:{λi}中不含特征值λi，其模為1且代數重數(algebraic multiplicity)大于幾何重數(geometric multiplicity)，即所有約當塊皆為一階.

以上判據可稱為經典譜判據.不僅當A非缺損時，而且當A缺損時，只要約當塊的階數有限，經典譜判據皆可用于判別強穩定性.經典譜判據的實用價值在于，對許多有意義的模型A為正規(normal)，此時使用這一判據的前提式(8)成立［5］.不過，對某些有意義的模型A為非正規，式(8)可能不成立并難以驗證［6］.這是經典譜判據的局限性.為消除這一局限性，Kreiss和Buchanon使用酉變換替換約當變換將任一方陣A展開為［13］

式中:U為J階酉陣，Γ為J階上三角矩陣:

式中:Γ的對角元素按模的大小作升序或降序排列.由此可得保證式(2)成立的充要譜判據

式中:K為與算子簇A無關的常數.式(13)提供了無前提條件的譜判據.上述譜判據可用于鑒別給定數值模擬方案的強穩定性.但是，有限模型的譜判據建立在排除所有內域節點控制方程與邊界條件聯合支持的失穩模態之上，因而不能揭示局部失穩的機理，即揭示某一節點與其鄰近節點相互作用產生的失穩模態.因此，對于出現了局部失穩的數值模擬方案，譜分析無助于為改進方案提供線索和論證局部穩定性.這是有限模型譜分析的局限性.

2.2 Godunov-Ryabenkii的工作和正則模態分析

Godunov和Ryabenkii首先通過簡單算例討論了有限模型譜分析引入的矛盾，然后提出消除矛盾的如下思想:將有限模型分解為若干基本的無限模型，并用正則模態分析研究用譜分析難以處理的無限模型的穩定問題［6］.這一工作開啟了穩定性研究的新階段.為了進一步揭示這項工作的意義，特別是對研究局部失穩問題的重要性，本節首先對文獻［6］所研究的算例作略為不同的分析，并介紹他們的和后續的主要研究結果;然后闡明此項成果為研究局部失穩問題所提供的啟示.

2.2.1 算例分析和相關研究結果

考察單向波動邊值問題:

就穩定性分析而言，式(15)的解可寫成式(1)的形式，其中，

A 的特征值為［14］

若將式(18)的解寫成式(1)的形式，其中，

依據2.1節提供的缺損矩陣A的譜判據可得Δτ*＜2，但依據強穩定性要求得到正確結果是Δτ*=1，再次出現矛盾.

為消除有限模型譜分析引入的矛盾，Godunov和Ryabenkii將上述有限模型分解為3個無限模型，即全無限和分別具有左、右端邊界條件的2個半無限模型，并提出保證有限模型穩定的如下G-R條件:不允許這3個無限模型中的任一個具有特征值的模大于1的模態［6］.由于在數值模擬方案的常規設計中通常已保證von Neumann穩定條件，問題歸結為保證半無限模型的穩定性.沿著這一方向后續研究工作表明G-R條件仍為穩定的必要條件，因為它未考慮當特征值的模為1時可能出現的失穩模態.為排除此模態需要補充一個條件，并針對一維一階雙曲型方程組證明了在一定附加條件下G-R條件加上補充的條件為穩定的充要條件，這個結果稱為GKS定理［7-9］.GKS定理所補充的條件的涵義如下:邊界條件和相鄰內域節點運動方程的離散格式皆不支持能量從邊界向內域流動的數值解，即離散格式不允許群速度從邊界指向內域［10］.

2.2.2 Godunov-Ryabenkii工作的評論

1)研究人員通過數值實驗早已觀察到局部失穩源于邊界與相鄰內域節點的相互作用，意識到需要將這一相互作用分離出來加以研究，并試圖實現這一分離以闡明失穩機理和制定消除局部失穩的措施［1-2，15］，但是，這些研究未能將這一“分離”的思想貫徹到底，未建立反映這一分離的基本數學模型，從而導致其研究結果的局限性.另一方面，在上述兩例中算子A的非正規性皆來自內域離散格式的非對稱性，即算例僅涉及內域穩定問題.實際上，若離散方案同時滿足von Neumann穩定條件，則有限模型的譜判據并不引入矛盾.后繼研究工作還表明，同時滿足有限模型譜判據和von Neumann穩定條件在許多情況下可以保證數值模擬方案的穩定性［15-16］.但是，這補充了附加條件的譜分析亦不能揭示邊界引入的局部失穩的機理.因此，Godunov等所研究的算例并未涉及本文提出的局部失穩問題.但是，他們隨后提出的用基本無限模型替換有限模型以及用正則模態分析闡明其失穩模態的思想對指導真實模型局部失穩機理的研究及尋求消除局部失穩的數值模擬方案具有普遍價值.這是因為，1)Godunov和Ryabenkii提出的基本無限模型完全實現了將控制局部失穩的因素從眾多影響因素中分離出來，使研究者的注意力集中到相互耦合的不多的幾類節點運動方程上.消除局部失穩的嘗試則歸結為改變這些方程或其中之一的結構或系數，以避免形成失穩模態的影響因素之間的匹配.這為探尋穩定方案提供了中肯的線索.2)他們首先在穩定性研究中使用的正則模態分析則為闡明真實模型的局部失穩機理提供了技術思路.

當k→0和h→0時，由上式知‖An‖趨于無窮.因此，在引入矛盾的這一Δτ區間內，經典譜判據不成立.其次，看第2個例子.將式(18)式給出的A展開成式(7)的形式，則

注意到 J隨 k→0而趨于無窮，對于在區間0＜Δτ＜1和1＜Δτ＜2內任一給定的 Δτ，式(8)皆不成立.更有意思的是，應用無前提條件的譜判據式(13)可得穩定準則0＜Δτ＜1.因此，譜判據并未引入矛盾.不過，瑕不掩瑜，這一評論僅為科學的非邏輯發展補充一例.

3 一維界面模型的穩定準則

為闡明正則模態分析的基本思想，本小節應用這一分析方法初步論證文獻［17］提出的一維分界面格式的經驗穩定準則，并作數值檢驗.考慮無限彈性桿中剪切波的傳播.設分界面位于x=0處，分界面兩側波動方程為

設分界面右側和左側步距分別為h和(c2/c1)h，離散 網 格 節 點 坐 標 為 xj=jh，j=0，1，2，…;xj=j(c2/c1)h，j=-1，-2，…(圖 1).分界面及其鄰近節點的遞推公式(文獻［14］式(15))，其二階格式為

圖1 一維非均勻網格示意圖Fig.1 One-dimensional non-uniform grid diagram

下面應用正則模態分析對格式(24)的穩定性進行初步論證，即該格式不支持特征值的模超過1的如下模態:

式中:φj(z)為非零模態向量，特征值z為復數.將式(26)代入式(24)并注意到式(25)給出的系數ξ1、η1、γ1、ζ1和 p=3/2 可得

式中:

差分方程式(27)的通解為

式中:κ為復數，φ為二維非零常數向量.將式(28)代入式(27)整理后得到

其中，式(27)成立的充要條件為

式(30)的展開式一般為κ的四階方程，但因四階項的系數和常數項為零而退化為二階方程.求解此二階方程得到 κ 的 2 個根 κ1、κ2，再將 κ =κ1，κ2代入式(29)得到對應的向量φ=φ1，φ2，其中，

由于內域穩定準則為0＜Δτ≤1，假定在此區間內 |z|＞ 1，不難證明|κ1|＜ 1，|κ2|＜ 1.注意到模態式(26)在界面節點兩側的表達式為

式中:b1、b2為常數.根據界面位移連續條件得到b1φ1(z)=b2φ2(z)，上式可寫成

因此，若格式(24)支持失穩模態式(26)，則將發生如式(32)所示失穩(圖2).

圖2 一維界面模型的失穩形式Fig.2 The form of G-R instability of the one-D model with an interface

關于式(24)支持上述失穩模態的論證是在|z|＞1的假定下做出的.但是，若采用前述經驗穩定準則，這一假定不成立.將式(32)代入遞推式(24)，若b1不為零得到

當p=3/2時從式(33)給出的2個相容的方程得到解z=-1.這與假定|z|＞1矛盾，即式(24)不支持|z|＞1的模態式(32).

設有限彈性桿中界面(x=0)右側和左側的波速分別為c1=100 m/s和c2=50 m/s，左端(x=-300 m)固定，右端(x=600 m)自由，界面阻抗比α=2.給定初始速度為零，初始位移f(x)為3次B樣條函數，其非零部分分布于區間［-25 m，0］.界面節點和內域節點的遞推公式由式(24)給出.兩端節點的遞推公式可作為界面節點遞推公式的極限情形得到，例如，若j=0的右側為自由面，令j=0和α=0，則式(24)退化為自由端遞推公式.若j=0的右側為固定端，令j=0和α→+∞，式(24)則退化為固定端遞推公式.為檢驗上述穩定準則，在數值實驗中僅變動界面節點和兩端點的p值，而對內域節點取用固定值p=3/2.實驗結果表明，當p≠3/2發生失穩.圖3給出界面節點的一組數值解(Δτ=1，界面左側和右側空間步距分別為3 m和6 m).

圖3 穩定方案與失穩方案的界面節點位移時程Fig.3 Displacement time history of interface node in stability plan and instability plan

4 弱穩定性分析

弱穩定性分析可用于線法(method of lines)建立的時空離散格式，亦可用于其他離散方法建立的數值解格式，只要當k→0時該格式退化為一組常微分方程，且在h→0時后者收斂于原偏微分方程問題的解.應用弱穩定性分析的主要問題是檢查常微分方程組式(5)初值問題的解是否收斂于連續模型的解.下面針對有限差分和有限元方法分別討論這一問題.

4.1 有限差分解法

就弱穩定性分析在有限差分解法中的應用而言，Strikwerda指出:若空間離散格式至少具有一階精度，則上述收斂性等價于如下適定條件［11］.

式中:C為與t無關的常數.假定式(5)中B為常系數方陣且維數不變，則適定條件式(34)可用B的譜判據檢驗［19］.但是，式(5)的收斂性涉及維數變化的一簇矩陣B，文獻［19］所提供的判據不再有效.除對十分簡單的模型可以證明式(34)成立，例如，Strikwerda通過類似Godunov和Ryabenkii采用的用無限模型代替有限模型的方法，對式(14)的一個簡單的空間離散格式證明了式(34)成立［11］.在一般情形下，式(34)是否成立則難以論證.因此，弱穩定性分析一般不能用于有限差分格式的穩定性分析.它的另一缺陷來自回避了實際的時空離散模型，從而排除了分析在實際數值計算過程中出現的某些現象的可能，例如，闡明局部失穩機理的可能.不過，在弱穩定性分析的前提成立的條件下，弱穩定性能消除所有可能的失穩形態，從而達到保證數值解格式穩定的目的.這是弱穩定性分析的價值所在，下面通過分析有限元解法闡明這一價值.

4.2 有限元解法

相對于具有局域近似特征的有限差分方法，基于泛函變分原理的有限元方法具有全局近似特征.這一特征使得用有限元方法建立的常微分方程組式(5)的收斂性易于得到保證.以線性粘彈性波的數值模擬為例，設邊界條件為給定位移、速度或應力，或給定應力與位移和速度的線性關系，則用有限元法得到的式(5)可寫成二階形式:

Mw..(t)+Cw·(t)+Kw(t)=0. (35)式中:M、C、K分別為質量、阻尼、剛度矩陣，其中，M對稱正定，K對稱半正定，C正定.Dupout證明了式(35)的解收斂于原連續模型的解［20］.除粘彈性波動外，淺水波等問題有限元離散所得常微分方程組的收斂性亦已獲得證明(參看文獻［20］引用的文獻).

對于若干具有實際意義的情形，式(35)時間離散格式的穩定性已得到證明.首先，若C為Rayleigh阻尼，則式(35)等價于解耦的一組單自由度方程，可采用常規譜判據判斷其穩定性［5］.其次，當C為一般正定陣時，用P層線性多步法離散式(35)，得

式中:ακ、βκ、γκ為實數，αp=1，Z 為時間平移算子，Zκvn=vn+κ.Gekeler將由 M、C、K 和 ρ(z)、σ(z)、τ(z)規定的離散格式的穩定問題歸結為［21］:是否存在實參數η∈-s，[]0，s＞0，使方程ρ(z)-ησ(z)=0的所有根z滿足|z|≤1且不存在|z|=1的復根?肯定回答是格式穩定的一個必要條件.在此條件下，穩定準則為

式中:λi為對稱陣M-1/2KM-1/2的特征值.

5 結束語

本文將Lax穩定性分析區分為強和弱2種穩定性分析方法.它們對進一步發展波動數值模擬技術各有其用途.強穩定性概念可用于任何時空離散模型的穩定性分析，特別適用于研究局部穩定性.它不僅可用于闡明實際模型中發生的典型局部失穩模態，而且對局部失穩模態取決于少數因素之匹配的認識，為探尋穩定的模擬方案提供了切實的技術路線.合理地使用弱穩定性分析可排除包含局部失穩在內的所有失穩模態，這一分析方法特別適用于完善波動的有限元模擬技術，例如，建立和論證穩定的、與內域數值解格式精度階匹配的實用人工邊界.

［1］LIAO Z P ，LIU J B.Numerical instabilities of a local transmitting boundary［J］.Earthquake Engineering ＆Structural Dynamics，1992，21(1):65-77.

［2］LIAO Z P，ZHOU Z H ，ZHANG Y H.Stable implementation of transmitting boundary in numerical simulation of wave motion［J］.Chinese Journal of Geophysics，2002，45(4):554-568.

［3］CRANK J，NICOLSON P.A practical method for numerical evaluation of solutions of partial differential equations of heat conduction type［J］.Proc Cambridge Philos Soc，1947，43:50-67.

［4］HOFFMAN J D.Numerical methods for engineers and Scientists［M］.New York:McGraw-Hill Inc，1992:657-691.

［5］HUGHES J R.Analysis of transient algorithms with particular reference to stability behavior［C］//Computational Methods for Transient Analysis.North-Holland:Elsevier Science，1983:67-155.

［6］GODUNOV S，RYABENKII V.Spectral stability criteria for boundary value problems for non-self-adjoint difference equations［J］.Russian Mathematical Surveys，1963，18(3):1-12.

［7］KREISS H O.Stability theory for difference approximations of mixed initial boundary value problems［J］.Mathematics of Computation，1968，22(104):703-714.

［8］GUSTAFSSON B，KREISS H O，SUNDSTR?M A.Stability theory of difference approximations for mixed initial boundary value problems［J］.Mathematics of Computation，1972，26(119):649-686.

［9］GUSTAFSSON B，KREISS H O，OLIGER J.Time dependent problems and difference methods［M］.New York:John Wiley＆ Sons，1995:496-620.

［10］TREFETHEN L.Group velocity interpretation of the stability theory of Gustafsson，Kreiss and Sundstr［J］.Journal of Computational Physics，1983，49(2):199-217.

［11］STRIKWERDA J C.Finite difference schemes and partial differential equations［M］.2nd ed.［S.1］:SIAM，2004:1-34.

［12］RICHTMYER R D，MORTON K W.Difference methods for initial-value problems［M］.2nd ed.New York:Interscience，1967:39-58.

［13］STRIKWERDA J C，WADE B A.A survey of the Kreiss matrix theorem for power bounded families of matrices and its extensions，in linear operators［C］//Linear Operators.Warszaw:Inst Math Pol Acad Sciences，1997:329-360.

［14］PARTER S V.Stability，convergence，and pseudo-stability of finite-difference equations for an over-determined problem［J］.Numerische Mathematik，1962，4(1):277-292.

［15］關慧敏，廖振鵬.局部人工邊界穩定性的一種分析方法［J］.力學學報，1996，28(3):376-380.

GUAN H M，LIAO Z P.A method for the stablitity analysis of local artificial boundaries［J］.Acta Mechanica Sinica，1996，28(3):376-380.

［16］SOUSA E.On the edge of stability analysis［J］.Applied Numerical Mathematics，2009，59(6)，1322-1336.

［17］謝志南，廖振鵬.波動方程數值模擬的一種顯式方法-界面節點遞推公式的構建［J］.力學學報，2011，43(1):154-161.

XIE Z N，LIAO Z P.An explicit method for numerical simulation of wave motion-constructing recursion formula for interface point［J］.Acta Mechanica Sinica，2011，43(1):154-161.

［18］廖振鵬，劉恒，謝志南.波動數值模擬的一種顯式方法——維波動［J］.力學學報，2009，41(3):350-360.

LIAO Z P，LIU H，XIE Z N.An explicit method for numerical simulation of wave motion:1-D wave motion［J］.Acta Mechanica Sinica，2009，41(3):350-360.

［19］HAIRER E，N?RSETT S P，WANNER G.Solving ordinary differential equations I:nonstiff problems［M］.2nd ed.Berlin:Springer，1993:80-89.

［20］DUPONT T.L2-estimates for Galerkin methods for second order hyperbolic equations［J］.SIAM Journal on Numerical Analysis，1972，10(5):880-889.

［21］GEKELER E.Linear multistep methods for stable differential equations［J］.Mathematics of Computation，1982，39(60):481-490.