王小花，潘文峰

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

王小花，潘文峰

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

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

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

0 引言

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

1 PML公式

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

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

類似于文獻［5，9－10］，設u^表示u的Laplace變換，定義如下:

引進坐標變換

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

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

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

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

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

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

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

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

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

接下來，引進一些輔助變量:

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

其中:

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

類似的推導可知，在二維空間中不存在變量ψ2，φ3，ζ3，此時PML方程簡化成如下形式:

其中:

2 穩(wěn)定性

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

考慮一個一般的偏微分方程的柯西問題:

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

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

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

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

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

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

其中:

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

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

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

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

其中:

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

3 數(shù)值實驗

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

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

對應的PML設置如下:

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

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

圖1 不同邊界條件的絕對誤差

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

4 結(jié)論

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

［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］楊彩虹，胡志興.一類具有飽和發(fā)生率的SEIR模型的穩(wěn)定性［J］.河南科技大學學報(自然科學版)，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］劉俊利，劉璐菊.具有媒體報道的傳染病模型穩(wěn)定性［J］.河南科技大學學報(自然科學版)，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

國家自然科學基金項目(11601402);湖北省自然科學基金項目(2014CFB865)

王小花(1990－)，女，河南信陽人，碩士生;潘文峰(1964－)，男，湖北天門人，教授，博士，碩士生導師，主要研究方向為偏微分方程數(shù)值解.

2016－12－21