時(shí)空多項(xiàng)式配點(diǎn)法求解三維Burgers方程*

曹艷華，張姊同，李 楠

（華東交通大學(xué) 理學(xué)院，南昌 330013）

引 言

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

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

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

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

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

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

1 時(shí)空多項(xiàng)式特解法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 配點(diǎn)圖：(a) 正方形區(qū)域規(guī)則取點(diǎn)；(b) 星形區(qū)域規(guī)則取點(diǎn)；(c) 正方形區(qū)域隨機(jī)取點(diǎn)；(d) 星形區(qū)域隨機(jī)取點(diǎn)(籃圈：內(nèi)點(diǎn)；紅星：邊界點(diǎn))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 時(shí)空區(qū)域中星形區(qū)域隨機(jī)取點(diǎn)（籃圈：內(nèi)點(diǎn)；紅星：邊界點(diǎn)）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í)驗(yàn)

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

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

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

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

本文的研究區(qū)域?yàn)檎襟w區(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ù)值實(shí)驗(yàn)中，我們將比較使用不同時(shí)刻 τ以及不同Reynolds數(shù)所得的數(shù)值模擬結(jié)果.如無特別說明，ni表示區(qū)域內(nèi)部的配點(diǎn)數(shù)，nd表示邊界??×(0,τ)∪?×{t=0} 上的配點(diǎn)數(shù)，nt表示測試點(diǎn)數(shù).

2.1 大τ 情形

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

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

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

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

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

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

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

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

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

圖6 τ=2時(shí)，正方體區(qū)域上系數(shù)矩陣的條件數(shù)：(a) 多項(xiàng)式特解系數(shù)矩陣的條件數(shù)；(b) 多項(xiàng)式基函數(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ù)值例子，可得下列結(jié)論.

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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