熱傳導(dǎo)方程的完美匹配層公式及其穩(wěn)定性分析

王小花，潘文峰

熱傳導(dǎo)方程的完美匹配層公式及其穩(wěn)定性分析

王小花，潘文峰

(武漢理工大學(xué)理學(xué)院，湖北武漢430070)

為了求解無(wú)界空間中的熱傳導(dǎo)方程，基于Laplace變換，引入若干個(gè)輔助變量，提出了一個(gè)熱傳導(dǎo)方程的完美匹配層(PML)公式。通過(guò)分析偏微分算子特征值實(shí)部的符號(hào)和特征向量的完備性，得到了PML方程的穩(wěn)定性。在二維空間中，常系數(shù)PML方程的柯西問(wèn)題是弱穩(wěn)定;在三維空間中，常系數(shù)PML方程的柯西問(wèn)題是強(qiáng)穩(wěn)定。數(shù)值實(shí)驗(yàn)結(jié)果表明:熱傳導(dǎo)方程PML公式的絕對(duì)誤差最大值大約是1.5×10－3，經(jīng)典Dirichlet邊界條件和Neumann邊界條件的絕對(duì)誤差最大值大約是2.5×10－2和3.0×10－2。因此，熱傳導(dǎo)方程PML公式可以顯著提高數(shù)值解的準(zhǔn)確性。

熱傳導(dǎo)方程;完美匹配層;穩(wěn)定性;輔助變量;Laplace變換

0 引言

熱量在無(wú)界介質(zhì)中的傳播是一個(gè)重要而又基本的問(wèn)題。對(duì)于波的傳播問(wèn)題，完美匹配層(perfectly matched layer，PML)算法［1］已經(jīng)被證明是一種有效且精確的方法。基于二維傅里葉變換，文獻(xiàn)［2］推導(dǎo)出了一個(gè)關(guān)于對(duì)流擴(kuò)散方程的PML公式，并且證明了在吸收層的外邊界，反射是非常小的。同時(shí)，數(shù)值實(shí)驗(yàn)也驗(yàn)證了完美匹配層算法的有效性和準(zhǔn)確性。對(duì)于波方程，目前有多種非常有效的PML算法，例如輔助微分方程完美匹配層(auxiliary differential equation perfectly matched layer，ADE-PML)算法［3］和復(fù)頻移完美匹配層(complex frequency shifted perfectly matched layer，CFS-PML)算法［4］。對(duì)于依賴(lài)時(shí)間的熱傳導(dǎo)方程，方程的穩(wěn)定性是非常重要的，不穩(wěn)定的數(shù)值方法可能會(huì)導(dǎo)致計(jì)算區(qū)域內(nèi)的數(shù)值解被污染。已經(jīng)有很多學(xué)者研究了偏微分方程的穩(wěn)定性［5－8］，然而，關(guān)于熱傳導(dǎo)問(wèn)題的完美匹配層算法的研究還比較少，特別是類(lèi)似于波方程的ADE-PML算法。本文基于Laplace變換，引入若干個(gè)輔助變量，嘗試推導(dǎo)出一個(gè)適用于熱傳導(dǎo)方程的ADE-PML算法，以期提高數(shù)值解的準(zhǔn)確性，并針對(duì)熱傳導(dǎo)方程的PML公式進(jìn)行穩(wěn)定性分析。

1 PML公式

考慮一個(gè)溫度場(chǎng)u在一個(gè)無(wú)界的三維空間中進(jìn)行傳導(dǎo)，并假設(shè)所有的熱源初始分布集中在一個(gè)給定的區(qū)域Ω=［－a1，a1］×［－a2，a2］×［－a3，a3］，a1，a2，a3＞0，因此，在區(qū)域Ω內(nèi)部，溫度場(chǎng)u(x1，x2，x3，t)滿(mǎn)足如下方程:

本文希望截?cái)酂o(wú)界區(qū)域并將計(jì)算限制在有限的區(qū)域Ω內(nèi)部。因此，假設(shè)Ω被包含于一個(gè)厚度為L(zhǎng)i(i=1，2，3)的PML，其中Li(i=1，2，3)為在xi(i=1，2，3)方向的厚度。在吸收層的內(nèi)部，溫度場(chǎng)u滿(mǎn)足一個(gè)修正的無(wú)熱源的以指數(shù)速度衰減的熱傳導(dǎo)方程。

類(lèi)似于文獻(xiàn)［5，9－10］，設(shè)u^表示u的Laplace變換，定義如下:

引進(jìn)坐標(biāo)變換

在區(qū)域Ω之外，u^滿(mǎn)足如下方程:

其中:ζi在吸收層xi＞ai是正的，在Ω內(nèi)部是0。若要求u^滿(mǎn)足如下的方程

顯然u在Ω的內(nèi)部保持不變，但是在吸收層內(nèi)呈指數(shù)衰減。

文獻(xiàn)［5］在Laplace變換域上給出了雙曲拋物系統(tǒng)標(biāo)準(zhǔn)的PML算法，但是涉及到了太多的輔助變量。本文嘗試使用較少的輔助變量，推導(dǎo)出一個(gè)新的關(guān)于熱傳導(dǎo)方程的PML公式。

由坐標(biāo)變換式(4)可知:

設(shè)γi=γi(ζi;s)，i=1，2，3，則:

根據(jù)式(6)可得:

根據(jù)式(7)，可以得到如下的恒等式:

結(jié)合式(8)和式(9)，可得:

接下來(lái)，引進(jìn)一些輔助變量:

使用Laplace逆變換將方程(10)變換回時(shí)域，可得:

其中:

在Ω內(nèi)部，吸收函數(shù)ζi(i=1，2，3)和輔助變量φ取值為0，因此，方程(12)又簡(jiǎn)化成了方程(1)?？梢钥吹?該P(yáng)ML公式只需要5個(gè)輔助變量φ1，φ2，φ3，ψ1，ψ2，并且沒(méi)有涉及到高階導(dǎo)數(shù)，算法的實(shí)現(xiàn)更加簡(jiǎn)單。

類(lèi)似的推導(dǎo)可知，在二維空間中不存在變量ψ2，φ3，ζ3，此時(shí)PML方程簡(jiǎn)化成如下形式:

其中:

2 穩(wěn)定性

本節(jié)中，將依次討論二維和三維情況下PML公式的穩(wěn)定性。對(duì)于拋物系統(tǒng)的穩(wěn)定性，文獻(xiàn)［11］有比較詳細(xì)的介紹，本文將引用文獻(xiàn)［11］中相關(guān)的理論證明方程(12)和方程(13)的穩(wěn)定性。

考慮一個(gè)一般的偏微分方程的柯西問(wèn)題:

根據(jù)文獻(xiàn)［11］可知:柯西問(wèn)題的適定性和穩(wěn)定性，可以通過(guò)算子P(iω)的特征值和特征向量來(lái)判斷。下面給出一個(gè)關(guān)于判別適定性和穩(wěn)定性的充要條件。

引理1［11］若存在不依賴(lài)于ω的常數(shù)C＞0(resp.C=0)使得算子P(iω)的所有特征值λ滿(mǎn)足:

則柯西問(wèn)題(14)是弱適定的(resp.弱穩(wěn)定的)。若對(duì)應(yīng)的特征向量是完備的，則柯西問(wèn)題(14)是適定的(resp.穩(wěn)定的)。

利用以上的穩(wěn)定性理論，可以得到下面兩個(gè)穩(wěn)定性定理。

定理1在二維空間中，對(duì)于任意的ζ1，ζ2＞0，若ζ1=ζ2，則PML公式的初值問(wèn)題是弱穩(wěn)定的。

證明設(shè)U=(u，φ1，φ2，ψ)T，可以將PML方程(13)寫(xiě)成如下形式:

其中:

設(shè)T=det(λI－P1(iω))，則:

由式(17)可知:只有ζ1=ζ2才能判斷PML方程的穩(wěn)定性。此外，當(dāng)ζ1=ζ2時(shí)，矩陣E不可對(duì)角化，因此PML方程是弱穩(wěn)定的。

定理2在三維空間中，對(duì)于任意的ζ1，ζ2，ζ3＞0，PML公式的初值問(wèn)題是強(qiáng)穩(wěn)定的。

證明令U=(u，φ1，φ2，φ3，ψ1，ψ2)，則可將式(12)寫(xiě)成如下形式:

其中:

因此，P(iω)的所有特征值的實(shí)部小于等于0，同時(shí)，矩陣E可對(duì)角化，故在三維空間中，對(duì)于任意的ζ1，ζ2，ζ3＞0，PML公式的初值問(wèn)題是強(qiáng)穩(wěn)定的。

3 數(shù)值實(shí)驗(yàn)

本節(jié)通過(guò)具體的數(shù)值算例說(shuō)明具有輔助微分方程的完美匹配層方法的有效性，并與經(jīng)典的狄利克雷(Dirichlet)邊界條件和紐曼(Neumann)邊界條件的結(jié)果相比較?？紤]如下的初值問(wèn)題:

其中:u0(x1，x2)=為初始的熱量分布，并且參數(shù)α=1。此時(shí)，初值問(wèn)題(20)具有解析解

對(duì)應(yīng)的PML設(shè)置如下:

其中:n=6;ζ=100;計(jì)算區(qū)域Ω=［－2，2］2;吸收層厚度L1=L2=0.3。對(duì)于對(duì)應(yīng)的PML方程，本文使用標(biāo)準(zhǔn)的有限元方法求解，采用線性基函數(shù)，網(wǎng)格的最大尺度是0.1，在PML吸收層的外邊界，采用齊次的Dirichlet邊界條件。

當(dāng)t=0.7時(shí)，不同的邊界條件所對(duì)應(yīng)的絕對(duì)誤差，如圖1所示。從圖1可以看出:PML邊界條件的絕對(duì)誤差最大值大約是1.5×10－3(見(jiàn)圖1a)，經(jīng)典的Dirichlet邊界條件和Neumann邊界條件的絕對(duì)誤差最大值分別約為2.5×10－2和3.0×10－2(見(jiàn)圖1b和圖1c)，因此，PML方法可以顯著提高數(shù)值解的準(zhǔn)確性。為了進(jìn)一步說(shuō)明PML邊界條件對(duì)原方程的影響，取t=0.7，x1=0，此時(shí)，數(shù)值解和解析解的差別如圖2所示。

圖1 不同邊界條件的絕對(duì)誤差

從圖2可以看出:在內(nèi)部區(qū)域，數(shù)值解和解析解幾乎完全重合;在吸收層內(nèi)部，u迅速衰減，到達(dá)吸收層的外部時(shí)幾乎為0，與波方程的結(jié)果類(lèi)似。

4 結(jié)論

本文推導(dǎo)出了一個(gè)熱傳導(dǎo)方程的PML公式，與文獻(xiàn)［5］中關(guān)于雙曲拋物系統(tǒng)的PML公式相比，本文的PML公式使用了較少的輔助變量。此外，證明了PML方程的柯西問(wèn)題的穩(wěn)定性。在二維空間中，常系數(shù)的完美匹配層方程的柯西問(wèn)題是弱穩(wěn)定的;在三維空間中，常系數(shù)的完美匹配層方程的柯西問(wèn)題是強(qiáng)穩(wěn)定的。通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證了熱傳導(dǎo)方程PML公式可以顯著提高數(shù)值解的準(zhǔn)確性。

［1］DRUSKIN V，GTTEL S，KNIZHNERMAN L.Near-optimal perfectly matched layers for indefinite Helmholtz problems［J］.Mathematics，2015，58(1):686－710.

［2］LANTOS N，NATAF F.Perfectly matched layers for the heat and advection-diffusion equations［J］.Journal of computational physics，2010，229(24):9042－9052.

［3］ATLURI S N，MARTIN R，KOMATITSCH D，et al.A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary dfferential equations(ADE-PML)［J］.Computer modeling in egineering＆sciences，2010，56(1):17－40.

［4］MA Y，CHEN S W，LU F，et al.The CFS-PML for 2-D WLP-FDTD method of dispersive materials［J］.International journal of applied electromagnetics and mechanics，2016，51(4):349－361.

［5］APPELD，HAGSTROM T，KREISS G.Perfectly matched layers for hyperbolic systems:general formulation，wellposedness，and stability［J］.SIAM journal on applied mathematics，2006，67(1):1－23.

［6］楊彩虹，胡志興.一類(lèi)具有飽和發(fā)生率的SEIR模型的穩(wěn)定性［J］.河南科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，2017，38(1): 78－83.

［7］COHEN G C，LIU Q H.Higher-order numerical methods for transient wave equations［J］.Springer，2002:349－354.

［8］劉俊利，劉璐菊.具有媒體報(bào)道的傳染病模型穩(wěn)定性［J］.河南科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，2016，37(2):88－91.［9］SJGREEN B，PETERSSON N A.Perfectly matched layers for Maxwell’s equations in second order formulation［J］.Journal of computational physics，2004，209(1):19－46.

［10］ABARBANEL S，GOTTLIEB D，HESTHAVEN J S.Long time behavior of the perfectly matched layer equations in computational electromagnetic［J］.Journal of scientific computing，2002，17(1):405－422.

［11］KREISS H O，LORENZ J.Initial-boundary value problems and the Navier-Stokes equations［M］.London:Academic Press

Inc，1989:44－62.

O193

A

1672－6871(2017)05－0070－06

10.15926/j.cnki.issn1672－6871.2017.05.015

國(guó)家自然科學(xué)基金項(xiàng)目(11601402);湖北省自然科學(xué)基金項(xiàng)目(2014CFB865)

王小花(1990－)，女，河南信陽(yáng)人，碩士生;潘文峰(1964－)，男，湖北天門(mén)人，教授，博士，碩士生導(dǎo)師，主要研究方向?yàn)槠⒎址匠虜?shù)值解.

2016－12－21