基于多重有限體積法解決含有耗散項(xiàng)的雙曲擴(kuò)散問(wèn)題?

張麗劍,羅躍生,高 洋

(1-哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，哈爾濱 150001;2-哈爾濱工程大學(xué)理學(xué)院，哈爾濱 150001)

1 引言

雙曲擴(kuò)散模型可看成是一種廣義的電報(bào)方程[1]，能夠模擬許多物理現(xiàn)象，例如在達(dá)西型多孔介質(zhì)中聲波的傳播以及平行流動(dòng)的粘性Maxwell流體的傳播[2]，模擬電壓和電流信號(hào)在同軸傳輸線中微小泄漏電導(dǎo)或電阻時(shí)的傳播[2]．當(dāng)然類(lèi)似的方程也同樣在不同領(lǐng)域中得到了應(yīng)用，例如模擬流體勢(shì)場(chǎng)中的擴(kuò)散過(guò)程[3]，含非線性阻尼項(xiàng)的運(yùn)輸記憶物理模型[4]和計(jì)算流體動(dòng)力學(xué)中的對(duì)流擴(kuò)散問(wèn)題的雙曲模型和各種熱傳導(dǎo)模型[5-7]，雖然在應(yīng)用中有很大的關(guān)聯(lián)，但方程引入對(duì)流項(xiàng)后，該問(wèn)題還未被徹底的研究．該問(wèn)題在某些初邊值條件下的解析解已經(jīng)給出，但至今無(wú)法得到一個(gè)廣義形式的解析解，甚至有些邊界為無(wú)界區(qū)域，因此研究其數(shù)值解是有必要的．

在本文中，考慮如下一維含耗散項(xiàng)的雙曲擴(kuò)散方程

其中P是常數(shù)，并且假設(shè)P的絕對(duì)值小于1，也就是說(shuō)P的絕對(duì)值小于擴(kuò)散項(xiàng)系數(shù)．因?yàn)椋魘P|大于擴(kuò)散項(xiàng)系數(shù)，從方程(1)的柯西問(wèn)題的漸進(jìn)性分析中能得到其解析解是不穩(wěn)定的[8]．

近些年，越來(lái)越多的中外學(xué)者研究了相對(duì)穩(wěn)定的偏微分方程數(shù)值算法和分析雙曲擴(kuò)散方程數(shù)值解方法[9-13]．本文針對(duì)問(wèn)題(1)，利用有限體積法思想，基于變限積分給出其構(gòu)造離散格式的方法；然后對(duì)其進(jìn)行了誤差估計(jì)和穩(wěn)定性分析．此外對(duì)其進(jìn)行了穩(wěn)定性分析，在第4節(jié)中，本文用能量分析法對(duì)離散格式進(jìn)行了先驗(yàn)估計(jì)和收斂性分析，得到該離散格式以二階收斂．最后針對(duì)不同初邊值條件給出數(shù)值算例，仿真結(jié)果表明數(shù)值結(jié)果跟理論分析能夠很好地吻合．

2 格式的構(gòu)造

對(duì)于問(wèn)題(1)，先將計(jì)算區(qū)域?=[a,b]進(jìn)行等步長(zhǎng)網(wǎng)格劃分，其中xi=ih,i=0,1,2,···,M,h=(b?a)/M；類(lèi)似地，tn=nτ,n=0,1,2,···,N,τ=T/N．記表示函數(shù)u在網(wǎng)格(ih,nτ)上的值．

首先，在控制區(qū)域?1=[xa,xb]×[ta,tb]內(nèi)，運(yùn)用有限體積法思想，對(duì)方程(1)兩邊作變限積分，其中xa,xb和ta,tb均分別為[a,b]和[0,T)中的變量

對(duì)(2)式化簡(jiǎn)可得

由于(3)式仍然存在偏導(dǎo)數(shù)的影響，因此繼續(xù)對(duì)上下限，在控制區(qū)域

內(nèi)，分別對(duì)ta,tb,xa和xb進(jìn)行積分

其中ε1,ε2,ε3和ε4分別為控制上下限積分區(qū)域大小的變限因子．

通過(guò)上述積分，將微分方程(1)完全轉(zhuǎn)化為積分方程(4)，接下來(lái)，本章采用不同點(diǎn)的拉格朗日插值對(duì)被積函數(shù)進(jìn)行插值逼近．選用二元九點(diǎn)拉格朗日插值法

對(duì)時(shí)間偏導(dǎo)項(xiàng)的u(x,t)進(jìn)行插值，其中Rl(x,t)為插值函數(shù)截?cái)嗾`差．而用二元六點(diǎn)拉格朗日插值法

即用時(shí)間n?1層和n+1層，空間i?1,i和i+1層數(shù)值對(duì)空間偏導(dǎo)項(xiàng)進(jìn)行插值，其中Rr(x,t)為插值函數(shù)截?cái)嗾`差．令εt=ε1=ε2,εx=ε3=ε4，可得

其中R1,R2,R2和R4分別為各項(xiàng)截?cái)嗾`差，具體表達(dá)式如下

由(18),(11)–(14)整理并舍去余項(xiàng)R1,R2,R2和R4，兩邊同乘以，即可得方程(1)的全離散格式，寫(xiě)成如下形式

顯然，全離散格式(15)需要兩層數(shù)值才能進(jìn)行運(yùn)算，在所給定的初始條件中，能夠得到第一層的精確值．對(duì)于第二層數(shù)值的逼近，本文通過(guò)對(duì)連續(xù)函數(shù)u(x,t)在t0處進(jìn)行Taylor展開(kāi)得到

此處

于是，能夠得到則(15)–(17)式，即為問(wèn)題(1)的全離散格式．

3 誤差估計(jì)

定理1假設(shè)u(x,t)足夠光滑，則全離散格式(15)的截?cái)嗾`差為

證明 在處理u時(shí)，針對(duì)等式(4)中，對(duì)時(shí)間偏導(dǎo)項(xiàng)，本文用了(xi?1,tn?1),(xi,tn?1),(xi+1,tn?1),(xi?1,tn),(xi,tn),(xi+1,tn),(xi?1,tn+1),(xi,tn+1)和(xi+1,tn+1)九個(gè)點(diǎn)做插值，所以由拉格朗日插值余項(xiàng)可得

其中ξ0,ξ1∈(xi?1,xi+1),δ0,δ1∈(tn?1,tn+1)．假設(shè)u(3,3)(ξ1,δ1),u(3,0)(ξ0,t),u(0,3)(x,δ0)有界，并取定ε1=ε2=C1τ,ε3=ε4=C2h，其中C1,C2∈(0,1)，而下文中出現(xiàn)的C為一般常數(shù)(即在不同地方表示不同的值)，則計(jì)算余項(xiàng)R1和R2，可得

對(duì)處進(jìn)行二元泰勒展開(kāi)，可得

由(20)可得

同理可得

此外

對(duì)處進(jìn)行泰勒展開(kāi)，可得

又因?yàn)棣?=ε2，則由(23)可得

結(jié)合(19),(21),(22)和(24)，可得

類(lèi)似地

對(duì)于空間偏導(dǎo)項(xiàng)，本文用了(xi?1,tn?1),(xi,tn?1),(xi+1,tn?1),(xi?1,tn+1),(xi,tn+1)和(xi+1,tn+1)六個(gè)點(diǎn)做插值，所以由拉格朗日插值余項(xiàng)可得

則由(27)計(jì)算可得

對(duì)處進(jìn)行二元泰勒展開(kāi)，可得

由(29)可得

同理可得

此外

對(duì)處進(jìn)行泰勒展開(kāi)，可得

又因?yàn)棣?=ε4，則由(32)可得

結(jié)合(28),(30),(31)和(33)，可得

類(lèi)似地

由于在對(duì)方程(1)進(jìn)行離散時(shí)，一共對(duì)x進(jìn)行了三次積分，對(duì)t進(jìn)行了三次積分，則等式兩邊應(yīng)除以h3τ3從而可得離散格式(15)的截?cái)嗾`差為

4 穩(wěn)定性分析

定理2當(dāng)時(shí)，全離散格式(15)是穩(wěn)定的．

證明 假設(shè)邊界條件精確滿(mǎn)足，采用Fourier分析法對(duì)全離散格式進(jìn)行穩(wěn)定性分析，用表示計(jì)算時(shí)產(chǎn)生的誤差，代入離散格式(15)，則滿(mǎn)足齊次方程

令σ是實(shí)數(shù)，η是復(fù)數(shù)，代入上式(36)，得到特征方程

其中

將η作等價(jià)替換，于是(37)可寫(xiě)成

由此可得，|η|<1的充要條件是Q1?Q2+Q3>0,Q1?Q3>0且Q1+Q2+Q3>0，對(duì)此，從特征方程(38)的系數(shù)中，令εt=C1τ,εx=C2h,C1,C2∈(0,1)，可得

因此，Q1?Q2+Q3>0，則

當(dāng)h→0時(shí)，Q1+Q2+Q3>0．由上述證明可知，離散格式(1)是無(wú)條件穩(wěn)定的．

5 離散格式的先驗(yàn)估計(jì)

方便起見(jiàn)，本文定義如下算子

則離散格式(7)–(11)忽略截?cái)嗾`差后，等式兩邊同乘以可改寫(xiě)成如下形式

整理(41)–(44)式可得

對(duì)于任意的本文定義離散空間內(nèi)積

定義離散L2空間范數(shù)

同樣定義離散Lp(1≤p≤∞)空間范數(shù)

當(dāng)p=∞，定義

引理1對(duì)任意有

證明 1) 因?yàn)?/p>

所以

2) 同理1)

3) 由2)可得

上述結(jié)論得證．

引理2(離散Gronwall不等式[14,15]) 設(shè)w(k),ρ(k)為非負(fù)網(wǎng)格函數(shù)，且ρ(k)為非遞減的．若

對(duì)所有的k成立，則w(k)6ρ(k)eCτk，對(duì)所有的k成立．

引理3(離散Sobelev不等式[14,15]) 對(duì)任意在有限區(qū)間[a,b]的離散范數(shù)0,1,2,···,M}，存在如下不等式

此處C(ε)和ε為常數(shù)，且C(ε)與ε有關(guān)．

定理3則離散格式(45)滿(mǎn)足如下等式

證明 將(45)與(un+1?un?1)作內(nèi)積，由邊界條件和引理1可得

其中

結(jié)合(48)–(54)式，整理可得

利用Young不等式

則(55)式可得

令

則(56)可寫(xiě)成

選擇足夠小的τ，使得1?Cτ=δ>0，則有

將(58)兩邊對(duì)0,1,2,···,N?1求和，可得

由引理2離散Gronwall不等式可得

由引理3離散Sobelev不等式可得

上述結(jié)論得證．

6 離散格式收斂性

定理4設(shè)則全離散格式(15)的解以∥·∥∞收斂到初邊值問(wèn)題(1)的解，且收斂階為

證明設(shè)為初邊值問(wèn)題(1)的解，且則由(45)式可得離散格式的截?cái)嗾`差Rn為

由定理1結(jié)論可對(duì)Rn進(jìn)行估計(jì)，設(shè)為離散格式(15)的解，定義由(59)可得

將(60)式與(en+1?en?1)作內(nèi)積，由邊界條件和引理1可得

令

則(61)可寫(xiě)成

選擇足夠小的τ，使得1?Cτ=δ>0，則有

將(62)兩邊對(duì)0,1,2,···,N?1求和，可得

由引理2離散Gronwall不等式可得

由此可得

上述結(jié)論得證．

7 數(shù)值算例

本節(jié)中，針對(duì)上述推導(dǎo)的離散格式進(jìn)行數(shù)值實(shí)驗(yàn)，通過(guò)比較精確解與數(shù)值解，來(lái)驗(yàn)證理論分析．令

此處u表示數(shù)值解，表示精確解，為了測(cè)量誤差和收斂速度，本章根據(jù)(46)和(47)定義∥e∥L2和∥e∥L∞．

考慮如下問(wèn)題[16,17]

初始條件

以及邊界條件

因?yàn)榭紤]到本問(wèn)題的物理域是一個(gè)無(wú)界區(qū)域，而用數(shù)值模擬時(shí)的計(jì)算域只能是一個(gè)有限區(qū)域，為實(shí)現(xiàn)數(shù)值目的，本節(jié)假設(shè)計(jì)算域是足夠大的，即在數(shù)值邊界處的解視為0為止(或者非常接近于0)，為了減小或避免數(shù)值邊界對(duì)計(jì)算結(jié)果的影響，本節(jié)給出的例子是T=50，即時(shí)間從0到40時(shí)刻的趨勢(shì)變化，選取的空間區(qū)域?yàn)閇?10,10]．

此處，令εt=τ/14,εx=h/14，取h=τ=0.1，計(jì)算時(shí)間域?yàn)閇0,50]，觀察當(dāng)P=0時(shí)，及無(wú)耗散項(xiàng)影響情況下的數(shù)值狀態(tài)，如圖1．

圖1: εt=τ/14,εx=h/14，取h=τ=0.1，當(dāng)T=20時(shí)數(shù)值圖及其俯視圖

圖2顯示的是，當(dāng)εt=τ/14,εx=h/14,h=τ=0.1，計(jì)算時(shí)間域?yàn)閇0,40]，P分別取?0.8,?0.5,?0.2,0.2,0.5,0.8時(shí)的數(shù)值俯視圖．從圖中可以觀察到在取不同的P時(shí)，該問(wèn)題解的速度場(chǎng)變化方向，并且當(dāng)|P|接近1時(shí)，即對(duì)流項(xiàng)系數(shù)越接近擴(kuò)散項(xiàng)系數(shù)時(shí)，波的傳播影響時(shí)間更長(zhǎng)，并且能夠影響接近邊界的值，解相對(duì)不穩(wěn)定．相反，當(dāng)|P|接近于0時(shí)，即該雙曲擴(kuò)散方程接近無(wú)耗散情況下，該波的傳播的能量能很快收斂，這跟文獻(xiàn)[8]分析的結(jié)果吻合．

圖2: P分別取?0.8,?0.5,?0.2,0.2,0.5,0.8時(shí)的數(shù)值俯視圖

8 結(jié)論

本文針對(duì)含耗散項(xiàng)的雙曲擴(kuò)散方程，利用基于變限積分的有限體積法，建立了一維的新的數(shù)值離散格式，詳細(xì)地介紹了構(gòu)造過(guò)程，并且進(jìn)行了誤差估計(jì)，其精度達(dá)到利用Fourier分析法對(duì)全離散格式進(jìn)行了穩(wěn)定性分析，發(fā)現(xiàn)所構(gòu)造的格式能在任意條件下都能達(dá)到穩(wěn)定．此外，對(duì)離散格式進(jìn)行了先驗(yàn)估計(jì)，討論了其收斂性，通過(guò)數(shù)值算例證實(shí)了上述結(jié)論．

參考文獻(xiàn)：

[1]陳宏霞，王同科.一維變系數(shù)對(duì)流擴(kuò)散方程第三邊值問(wèn)題的緊有限體積方法[J].工程數(shù)學(xué)學(xué)報(bào)，2014,31(6):889-902 Chen H X,Wang T K.Compact finite volume method for 1D convection dif f usion equations with variable coefficients and third boundary conditions[J].Chinese Journal of Engineering Mathematics,2014,31(6):889-902

[2]荊菲菲，蘇劍，陳浩.二階橢圓邊值問(wèn)題的一種新型穩(wěn)定化非協(xié)調(diào)混合有限元方法[J].工程數(shù)學(xué)學(xué)報(bào)，2013,30(6):846-854 Jing F F,Su J,Chen H.A new stabilized nonconforming-mixed finite element method for the second order elliptic boundary value problem[J].Chinese Journal of Engineering Mathematics,2013,30(6):846-854

[3]王強(qiáng)，陶建華.非線性擬雙曲方程的有限元配置法數(shù)值分析[J].工程數(shù)學(xué)學(xué)報(bào)，2012,29(1):96-106 Wang Q,Tao J H.Numerical analysis of the finite element collocation method for the nonlinear pseudohyperbolic equation[J].Chinese Journal of Engineering Mathematics,2012,29(1):96-106

[4]Mohanty R K.New unconditionally stable dif f erence schemes for the solution of multi-dimensional telegraphic equations[J].International Journal of Computer Mathematics,2009,86(12):2061-2071

[5]Biazar J,Eslami M.Analytic solution for telegraph equation by dif f erential transform method[J].Physics Letters A,2010,374(29):2904-2906

[6]Ara′ujo A,Das A K,Neves C,et al.Numerical solution for a non-Fickian dif f usion in a periodic potential[J].Communications in Computational Physics,2013,13(2):502-525

[7]Mac′?as-D′?az J E.Sufficient conditions for the preservation of the boundedness in a numerical method for a physical model with transport memory and nonlinear damping[J].Computer Physics Communications,2011,182(12):2471-2478

[8]G′omez H,Colominas I,Navarrina F,et al.A discontinuous Galerkin method for a hyperbolic model for convection-dif f usion problems in CFD[J].International Journal for Numerical Methods in Engineering,2007,71(11):1342-1364

[9]Kulish V,Poletkin K V.A generalized relation between the local values of temperature and the corresponding heat flux in a one-dimensional semi-infinite domain with the moving boundary[J].International Journal of Heat and Mass Transfer,2012,55(23-24):6595-6599

[10]Lin H C,Char M I,Chang W J.Soret ef f ects on non-Fourier heat and non-Fickian mass dif f usion transfer in a slab[J].Numerical Heat Transfer,Part A:Applications:An International Journal of Computation and Methodology,2009,55(12):1096-1115

[11]Zauderer E.Partial Dif f erential Equations of Applied Mathematics[M].New York:John Wiley&Sons Inc.,1989

[12]Zhang Z Y,Deng D W.A new alternating-direction finite element method for hyperbolic equation[J].Numerical Methods for Partial Dif f erential Equations,2007,23(6):1530-1559

[13]Fernandes R,Fairweather G.A Crank-Nicolson and ADI Galerkin method with quadrature for hyperbolic problems[J].SIAM Journal on Numerical Analysis,1991,28(7):1265-1281

[14]Zhang Z Y,Deng D W.A new unconditionally stable ADI compact scheme for the two-space-dimensional linear hyperbolic equation[J].International Journal of Computer Mathematics,2010,87(10):2259-2267

[15]Wang T C,Guo B L,Zhang L M.New conservative dif f erence schemes for a coupled nonlinear Schr¨oinger system[J].Applied Mathematics and Computation,2000,217(4):1604-1619

[16]Fernandes R I,Fairweather G.An alternating direction Galerkin method for a class of second-order hyperbolic equations in two space variables[J].SIAM Journal on Numerical Analysis,1991,28(5):1265-1281

[17]Ara′ujo A,Neves C,Sousa E.An alternating direction implicit method for a second-order hyperbolic dif f usion equation with convection[J].Applied Mathematics and Computation,2014,239(2):17-28