曹艷華，張姊同，李 楠

（華東交通大學 理學院，南昌 330013）

引 言

Burgers 方程可以用來模擬沖擊波的傳播和反射，是一個被廣泛應用于流體力學、非線性聲學、氣體動力學等領域的非線性偏微分方程[1-3].由于Burgers 方程是從Navier-Stokes 方程簡化而來，因此對于該方程的研究眾多，其中不乏非常有意義的研究成果.Hopf[2]和Cole[4]提出了著名的Hopf-Cole 變換來研究Burgers 方程的基本性質(zhì).他們發(fā)現(xiàn)Reynolds 數(shù)對于平衡方程中的非線性對流項和黏性擴散項起著舉足輕重的地位.當Reynolds 數(shù)變大的時候非線性對流項占主導地位，由此導致了Burgers 方程的強非線性性，使得其求解異常困難.盡管通過Hopf-Cole 變換可以求得一些特殊類型Burgers 方程的精確解，但對于絕大多數(shù)Burgers 方程而言，其精確解是難以得到的.學生在學習Burgers 方程的強非線性理論時，不易理解其含義[5-7].

Gulsu 使用有限差分法求解具有小Reynolds 數(shù)的Burgers 方程[8]，對于具有大Reynolds 數(shù)的Burgers 方程，有限差分法的求解精度很差[9].Turgut 等采用半步長方法[10]，即時間方向采用半隱式有限差分格式、空間方向采用漸進展開格式求解具有大Reynolds 數(shù)的Burgers 方程.除了采用傳統(tǒng)的網(wǎng)格類方法，近些年來，無網(wǎng)格類方法也被廣泛用來求解Burgers 方程[11-17].

所謂時空類方法，即在求解發(fā)展方程時將時間變量t視為普通的空間變量x，其中著名的時空類方法有時空有限元法[18]、間斷Galerkin 有限元法[19]以及Mohammed 等提出的局部時空徑向基函數(shù)等方法[20].

2017年，Dangal 等提出了一種新的多項式特解法用于求解偏微分方程[21].2020年，Cao 等在多項式特解方法的基礎上，提出了一種時空多項式特解法求解一維和二維的發(fā)展方程[22-23].受其啟發(fā)，本文采用時空多項式特解法求解三維具有大Reynolds 數(shù)的Burgers 方程.相比于其他的無網(wǎng)格類方法，多項式特解配點法只需知道配點的位置信息，別的任何未知參數(shù)都不需要確定，其簡便性遠遠超過現(xiàn)有的其他時空類無網(wǎng)格方法，例如，使用時空徑向基函數(shù)方法時，不僅需要確定徑向基函數(shù)的形狀參數(shù)，還需要確定配點的中心位置、個數(shù)以及插值半徑，這些參數(shù)的確定沒有一般規(guī)律可循.數(shù)值例子表明，隨著Reynolds 數(shù)的增加，時空多項式特解法不僅可以保持數(shù)值解的高精確性，其穩(wěn)定性也沒有受到任何影響.

注1所有以多項式作為基函數(shù)的方法中，由于系數(shù)矩陣隨著多項式階的增加，其條件數(shù)迅速增大，很快數(shù)值溢出，因此本文的所有算例均采用了多尺度技巧以降低系數(shù)矩陣的條件數(shù)，得到穩(wěn)定的數(shù)值解[22].

注2多尺度技巧可以有效降低多項式特解中系數(shù)矩陣的條件數(shù)，但對于多項式直接作為基函數(shù)得到的系數(shù)矩陣，多尺度技巧不能有效降低系數(shù)矩陣的條件數(shù)，這也是多項式函數(shù)直接作為基函數(shù)求解三維Burgers 方程得到的近似解精度不高的原因之一.

1 時空多項式特解法

考慮如下的三維偏微分方程初邊值問題：

其中，Ω為三維空間中一有界連通的規(guī)則或不規(guī)則區(qū)域，f(x,t),g(x,t)和h0(x)為給定函數(shù)，?t=?/?t，L為線性或非線性偏微分算子，B為邊界算子，(0,τ]為待求時間區(qū)間.

在時空類方法中，需要將三維待求解問題轉化為相應的四維問題，新的“空間變量”為“t+(x,y,z)”，將?×{t=0} 和 ??×(0,τ)作為新四維問題的邊界.在時空類方法中，不再區(qū)分時間變量t和 空間變量x=(x,y,z)，因此，相應于方程(1)的新四維問題為

在新四維空間中，d階完全多項式基函數(shù)為

其中 Λ={0≤i+j+k+l≤d}.

空間 Qd中，線性無關的基函數(shù)個數(shù)為(d+1)(d+2)(d+3)(d+4)/24.假設方程(2)的解h(x,t)可以寫為特解 χijkl(x,t)的一個線性組合：

其中 ζijkl為待求解系數(shù)，特解 χijkl(x,t)滿足

L1為L 中的空間常系數(shù)線性算子部分.在四維空間中，關于二階線性常系數(shù)偏微分方程的特解 χijkl(x,t)有以下引理.

引理1[21]令(ω1,ω2,ω3,ω4)=(x,y,z,t)，則對于下列偏微分方程

為式(6)的一個多項式特解，其中d=i+j+k+l,L=

在區(qū)域 ? 內(nèi)部任取ni個配點，邊界 ? ?×(0,τ)上任取nd1個配點，邊界? ×{t=0} 上任取nd2個配點，配點總數(shù)目為N=ni+nd1+nd2，則方程(2)的配點型強格式為

迭代求解代數(shù)方程組(8)，可得到待定系數(shù){ζijkl|Λ}.將{ζijkl|Λ}代入式(4)即得到方程(2)的近似解.

為了比較說明多項式特解法的優(yōu)點，假設方程(2)的解h(x,t) 可表示為多項式基函數(shù)?ijkl(x,t)=xiyjzktl的一個線性組合：

其中 ?ijkl為待定系數(shù).將(x,t) 代入方程(2)，且在區(qū)域 ? 內(nèi)部任取ni個配點，邊界 ? ?×(0,τ)上任取nd1個配點，邊界?×{t=0}上任取nd2個配點，配點總數(shù)目為N=ni+nd1+nd2，則方程(2)的配點型強格式為

迭代求解代數(shù)方程組(10)，得到待定系數(shù){?ijkl|Λ}.將{?ijkl|Λ}代入式(9)即得到方程(2)的近似解.

注3為使代數(shù)方程組(8)和(10)可解，配點的數(shù)目N 需大于完全多項式空間 Qd的維數(shù)(d+1)(d+2)(d+3)(d+4)/24.

本文的計算區(qū)域為規(guī)則區(qū)域和不規(guī)則區(qū)域，所用配點有兩種方式：一是內(nèi)點和邊界點均為規(guī)則配點，二是內(nèi)點和邊界點均為隨機配點.以二維區(qū)域為例，圖1(a)、1(b)所示的是規(guī)則配點，圖 1(c)、1(d)所示的是隨機配點.二維非規(guī)則空間區(qū)域在新的三維時空域中隨機取點如圖2所示，規(guī)則二維空間區(qū)域在新的三維時空域中隨機取點類似.

圖1 配點圖：(a) 正方形區(qū)域規(guī)則取點；(b) 星形區(qū)域規(guī)則取點；(c) 正方形區(qū)域隨機取點；(d) 星形區(qū)域隨機取點(籃圈：內(nèi)點；紅星：邊界點)Fig.1 The collocation point diagram: (a) the regular domain with regular points; (b) the irregular domain with regular points; (c) the regular domain with random points; (d) the irregular domain with random points (blue circles: interior collocation points; red stars: boundary collocation points)

圖2 時空區(qū)域中星形區(qū)域隨機取點（籃圈：內(nèi)點；紅星：邊界點）Fig.2 An irregular domain with random points in the space-time domain (blue circles: interior collocation points; red stars: boundary collocation points)

2 數(shù)值實驗

數(shù)值解與精確解之間的平方根誤差定義為

考慮下列三維Burgers 方程的初邊值問題：

其中，μ=1/Re＞0為黏性系數(shù)，Re為Reynolds 數(shù).方程(11)的精確解為

f(x,t),g(x,t)和h0(x)為基于精確解(12)的函數(shù).

本文的研究區(qū)域為正方體區(qū)域和bumpy-shaped 型非規(guī)則區(qū)域，如圖3所示.

圖3 空間求解區(qū)域：(a) 正方體區(qū)域；(b) bumpy-shaped 區(qū)域Fig.3 The spatial solution domains: (a) the cubic domain; (b) the bumpy-shaped domain

在數(shù)值實驗中，我們將比較使用不同時刻 τ以及不同Reynolds數(shù)所得的數(shù)值模擬結果.如無特別說明，ni表示區(qū)域內(nèi)部的配點數(shù)，nd表示邊界??×(0,τ)∪?×{t=0} 上的配點數(shù)，nt表示測試點數(shù).

2.1 大τ 情形

本小節(jié)的算例采用時空多項式配點法和直接使用多項式基函數(shù)時空配點法求解Reynolds 數(shù)Re=1時的方程(11).

表1給出了 τ=1，配點ni=3 819,nd=5 220,nt=1 919時，兩種區(qū)域?qū)Σ煌A多項式特解的計算時間和平方根誤差.可以看出，隨著特解階的增加，平方根誤差迅速下降；特解的階越高，計算所需時間越長.同時，可以看出，立方體區(qū)域上的平方根誤差下降得比bumpy 型區(qū)域上的平方根誤差快.

表1 ni=3 819,nd=5 220,nt=1 919,τ=1時，兩種區(qū)域多項式特解的平方根誤差Table 1 The RMSE with ni=3 819,nd=5 220,nt=1 919,τ=1 for cubic and bumpy-shaped domains

表2給出了 τ=1，配點ni=3 819,nd=5 220,nt=1 919時，兩種區(qū)域?qū)Σ煌A多項式基函數(shù)的計算時間和平方根誤差.可以看出，隨著多項式階的增加，平方根誤差未有明顯下降；多項式基函數(shù)的階越高，計算所需時間越長.通過分析多項式基函數(shù)的系數(shù)矩陣，可以看到多尺度方法對其系數(shù)矩陣的條件數(shù)未有明顯改進，這也是多項式基函數(shù)的平方根誤差并不隨著多項式基函數(shù)階的增加而降低的原因之一.

表2 ni=3 819,nd=5 220,nt=1 919,τ=1時，兩種區(qū)域多項式基函數(shù)的平方根誤差Table 2 The RMSE with ni=3 819,nd=5 220,nt=1 919,τ=1 for cubic and bumpy-shaped domains

τ=2的數(shù)值結果如圖4所示.易見，對這兩種區(qū)域，當特解基函數(shù)的階小于或等于12 時，所得計算結果的誤差都以指數(shù)量級迅速下降.此后，再增加特解基函數(shù)的階，所得計算結果的精度不再有顯著提高，但也不會降低.當 τ=100時，兩種區(qū)域的計算結果如圖5所示.顯然，當特解基函數(shù)的階小于或等于12 時，所得計算結果都快速收斂.

圖4 τ=2,Re=1時，兩種區(qū)域的計算結果Fig.4 Results for τ=2 andRe=1

圖5 τ=100,Re=1時，兩種區(qū)域的計算結果Fig.5 Results for τ=100 andRe=1

為進一步分析所得結果，τ=2時，正方體求解區(qū)域上的系數(shù)矩陣的條件數(shù)如圖 6(a)所示.顯然，在多項式特解中，使用多尺度技巧與否，所得矩陣的條件數(shù)相差巨大.如果不使用多尺度技巧，其條件數(shù)隨著多項式特解階的增加迅速增大，很快便數(shù)值溢出：多項式特解的階等于12 時，所得數(shù)值結果發(fā)散.多項式特解的階等于20 時，條件數(shù)溢出.使用多尺度技巧，可以有效控制條件數(shù)的增長，獲得穩(wěn)定的解.因此，在多項式類的解法中，多尺度技巧的使用是至關重要的.Bumpy-shaped 區(qū)域上條件數(shù)的情形與正方體求解區(qū)域上的情形類似.對于多項式基函數(shù)而言，系數(shù)矩陣的條件數(shù)使用和不使用多尺度技巧差別不大，如圖 6(b)所示，這或許是直接使用多項式基函數(shù)所得的近似解精度不高的原因之一.

圖6 τ=2時，正方體區(qū)域上系數(shù)矩陣的條件數(shù)：(a) 多項式特解系數(shù)矩陣的條件數(shù)；(b) 多項式基函數(shù)系數(shù)矩陣的條件數(shù)Fig.6 Matrix condition numbers of the cubic domain for τ=2: (a) the condition number of polynomial particular solutions; (b) the condition number of polynomial basis functions

根據(jù)本小節(jié)數(shù)值例子，可得下列結論.

1） 對于規(guī)則或不規(guī)則的計算區(qū)域，采用不同數(shù)目的配點，計算結果會有所不同，但是這些結果之間的差別很小，不足以影響到計算精度.

2） 采用不同數(shù)目配點的最大的區(qū)別在于計算時間，配點越多，所需的計算時間也越長.以正方體區(qū)域為例，ni=6 500，nd=6 000，nt=256時 所需的計算時間約為8 930 s.如果ni增加到10 000，nd增加到8 128，nt保持不變，則所需的計算時間增加4 529 s.由此可見，配點數(shù)目嚴重影響計算時間.在保證近似解精度的情況下，應選用適當數(shù)目的配點.對于非規(guī)則區(qū)域，其情形與正方體區(qū)域的情形類似，此處不再贅述.

3） 本算法中的多項式特解基函數(shù) χijkl需要提前生成.生成1 階到30 階的多項式特解，所需總時間約為10 415 s.值得慶幸的是，特解基函數(shù) χijkl一旦生成便可重復使用，甚至對于不同的微分方程，只要其中的線性微分算子項相同，也可使用.

2.2 不同Reynolds 數(shù)的情形

本小節(jié)中所有算例的τ均取1.在求解Burgers 方程時，的取值對所求數(shù)值解的精度影響巨大.現(xiàn)有研究結果表明，Reynolds 數(shù)Re越大，所得數(shù)值解的精度越低，求解難度越大.

當 μ=1，所得結果如圖7所示.正方體區(qū)域上數(shù)值解的精度達到10-11，即使在非規(guī)則區(qū)域bumpy-shaped上，數(shù)值解的精度仍然能達到10-9.當μ=0.5時，所得結果如圖8所示.正方體區(qū)域上數(shù)值解的精度達到10-8，非規(guī)則區(qū)域bumpy-shaped 上數(shù)值解的精度約為1 0-6.繼續(xù)增大Reynolds 數(shù)，μ=0.1，正方體區(qū)域上數(shù)值解的精度達到10-5，而非規(guī)則區(qū)域bumpy-shaped 上數(shù)值解的精度約為1 0-3，且數(shù)值解不再像 μ=1時那樣光滑.因此，本文數(shù)值結果進一步驗證了Reynolds 數(shù)越大，所得數(shù)值解的精度越低，求解難度越大，這與以往的研究結果一致.

圖7 τ=1,μ=1時，兩種區(qū)域的計算結果Fig.7 Results for τ=1 andμ=1

圖8 τ=1,μ=0.5時，兩種區(qū)域的計算結果Fig.8 Results for τ=1 andμ=0.5

表3給出了 τ=1，配點ni=3 819,nd=5 220,nt=1 919,μ=0.5時，兩種區(qū)域直接使用多項式基函數(shù)的計算時間和平方根誤差.可以看出，隨著多項式基函數(shù)階的增加，平方根誤差幾乎沒有降低.

表3 ni=3 819,nd=5 220,nt=1 919,μ=0.5時，兩種區(qū)域多項式基函數(shù)的平方根誤差和計算所需時間Table 3 The RMSE with ni=3 819,nd=5 220,nt=1 919,μ=0.5 for cubic and bumpy-shaped domains

3 結 論

三維Burgers 方程因其維數(shù)高，大Reynolds 數(shù)時的強非線性性，現(xiàn)有研究結果不多.本文提出了一種新型時空多項式特解配點法求解三維Burgers 方程，不同于以往的多項式函數(shù)直接做基函數(shù)的方法，多項式特解做基函數(shù)的優(yōu)點是對控制方程中的線性導數(shù)項不進行求導，因此多項式特解的階比較高時，所得系數(shù)矩陣中元素的量級差別不大，從而算法的穩(wěn)定性強，精度很高.采用配點法時，時空方法的維數(shù)比初始問題的維數(shù)增加了一維，所需配點數(shù)目較多，下一步的工作可以通過采取有效的降階技巧減少對配點的要求，將系數(shù)矩陣變?yōu)橄∈杈仃?，以期得到更好的模擬結果.通過本文的工作，使學生對高維強非線性的Burgers 方程的數(shù)值求解有更加深入的理解.