Darcy-Brinkman方程的殘量型后驗誤差估計?

張秋雨，黃鵬展

(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆烏魯木齊830046)

0 引言

在一些物理和工程的實際數(shù)值模擬過程中,我們經(jīng)常會遇到局部奇異、邊界層和近似斷裂問題.當(dāng)我們用有限元方法求解這類問題時,這些局部奇異的地方往往會在很大程度上影響全局的計算精度.自適應(yīng)有限元方法(簡稱AFEM)[1–3]是解決以上問題的一種比較常用的方法,主要思想是細(xì)化局部網(wǎng)格或者提升基函數(shù)的自由度.對于有限元方法來說,進(jìn)行網(wǎng)格局部加密是比較容易實現(xiàn)的,但如何確定這些關(guān)鍵的局部區(qū)域是我們需要重點關(guān)注的.后驗誤差估計一般分為基于方程殘量型的后驗誤差、基于梯度重構(gòu)型的后驗誤差和基于投影算子型的后驗誤差,根據(jù)這些誤差我們可以構(gòu)造對應(yīng)的誤差估計子,通過對這些估計子的計算和處理,我們可以用其指示出原網(wǎng)格需要加密的區(qū)域,通過加密策略進(jìn)而得到加密后的網(wǎng)格.

本文主要研究Darcy-Brinkman方程的穩(wěn)定化有限元方法的殘量型后驗誤差估計.其中ν代表粘性系數(shù),??=ΓT+ΓB是多邊形區(qū)域邊界,u,p,T,C分別代表速度場、壓力場、溫度場和濃度場,Da表示Darcy數(shù),Dc表示質(zhì)量擴(kuò)散系數(shù)，g是重力加速度,βT和βC分別為熱擴(kuò)散系數(shù)和物質(zhì)的溶解系數(shù),γ表示熱擴(kuò)散率.關(guān)于Darcy-Brinkman方程的先驗誤差估計的理論可以參考[4-6].

1 基于穩(wěn)定化有限元方法的殘量型后驗誤差估計子

本節(jié)主要介紹穩(wěn)定化有限元方法以及殘量型后驗誤差估計子的一些基本概念、推導(dǎo)過程,詳細(xì)內(nèi)容可以參見文獻(xiàn)[7-10].

考慮如下函數(shù)空間

定義Banach 空間對應(yīng)的范數(shù)這樣方程(1)的弱形式可寫為求解使得

當(dāng)采用有限元進(jìn)行離散這個問題時,我們考慮最低階的協(xié)調(diào)有限元對P1/P0/P1/P1(線性速度/常數(shù)溫度/線性溫度/線性濃度),選取如下的有限元子空間

這里的τh={K}是擬一致的網(wǎng)格剖分.為了使離散問題穩(wěn)定、滿足離散的inf-sup條件,我們需要采取穩(wěn)定化方法.在這里我們選用基于壓力投影的穩(wěn)定化方法[11],則離散變分問題可以寫成

穩(wěn)定項S(ph,qh)=(ph?Πph,qh?Πqh).Π定義為:L2(?)→P1(K),可參見文獻(xiàn)[11].

定義殘量泛函R:Mh×Qh×Wh×Lh→R,

根據(jù)殘量泛函,可以看出殘量包含單元殘量和邊界殘量兩部分,單元殘量反應(yīng)了有限元解在網(wǎng)格單元內(nèi)部逼近原始方程的程度,而邊界殘量在這里即單元間殘量,反映了法向通量在網(wǎng)格單元邊界的跳躍.定義∈h為所有邊的集合,K為τh的元素,hK為單元K的最長邊,he為單元K的邊界邊,ne為邊上的法向量.由以上的殘量泛函來構(gòu)建殘量型估計子

則全局估計子為

2 殘量型后驗誤差估計子的誤差估計

為了能應(yīng)用Ver¨ufrth和Ainsworth對非線性偏微分方程建立起的后驗誤差估計的理論框架[7–9],我們定義C?(X,Y)是Banach 空間X到Y(jié)的連續(xù)線性映射,并裝備范數(shù)‖·‖C?(X,Y).定義Isom(X,Y)?C?(X,Y)是X到Y(jié)的線性同胚的一個子集.XD?X,Y?D?Y?.Y?是Y的對偶空間,Y?=C?(Y,R).F∈C1(X,Y?)是一個連續(xù)的可微函數(shù),Fh滿足

DF(U?)是F對U?的線性化算子.

定理1[7,9,12]記Uh∈Xh是方程Fh(Uh)=0的數(shù)值解,且‖F(xiàn)h(Uh)‖Y?

h無限小.假定存在一個限制算子Rh∈C?(Y,Yh),Y→Y的恒等算子IdY,Y的一個有限維子空間?Yh=span{ψi},且?Fh:Xh→?Y?h是F在Uh的近似.則有

定理2記[u?,p?,T?,C?]是弱形式(2)的解,[uh,ph,Th,Ch]∈[Mh,Qh,Wh,Lh]是離散形式(3)的解且充分接近[u?,p?,T?,C?],并滿足定理1的條件,則我們有

證明首先,我們要驗證F的可微性,并在[u?,p?,T?,C?]的鄰域內(nèi)是Lipschitz連續(xù).令DF?∈C?(X,Y?),

我們可以看出F在[u?,p?,T?,C?]是可微的.

令DF1(·)是[u1,p1,T1,C1]處的導(dǎo)數(shù),對任意的[v,q,t,r]∈Y,[u,p,T,C]∈B([u?,p?,T?,C?],R?),有

則

即F是Lispschitz連續(xù)的.

根據(jù)和F的定義,我們可以得到

又因為

我們有

由(4)-(9),我們可以得到

3 數(shù)值實驗

在給出數(shù)值算例之前,我們需要引入一些記號.iter是自適應(yīng)加密的迭代次數(shù);Mesh是網(wǎng)格的單元個數(shù);uH1是速度數(shù)值解的H1相對誤差;pL2是壓力數(shù)值解的L2相對誤差;TH1是溫度數(shù)值解的H1相對誤差;CH1是濃度數(shù)值解的H1相對誤差.

3.1 真解算例

考慮方形區(qū)域?=[0,1]×[0,1],我們給出以下真解,其中r1=4.1,r2=4.1.

首先我們先給出一致網(wǎng)格的有限元數(shù)值模擬結(jié)果,來作為基于以上估計子的自適應(yīng)算法的對比.從表1和表2的最后一次加密結(jié)果我們可以看出,用自適應(yīng)加密策略生成的網(wǎng)格來進(jìn)行計算時比一致網(wǎng)格的效率更高.在表2中的第五次迭代時,我們只用了7 511個單元,就得到比一致網(wǎng)格8 192個單元低一倍的相對誤差,體現(xiàn)出了新方法的有效性.另外圖1描繪出了相關(guān)的網(wǎng)格圖和數(shù)值結(jié)果.

表1 一致網(wǎng)格的數(shù)值誤差

表2 自適應(yīng)網(wǎng)格的數(shù)值誤差

圖1 第一行是自適應(yīng)加密的得到網(wǎng)格,第二行在自適應(yīng)網(wǎng)格下得到的數(shù)值解uh,ph,Th

圖2 第一行是自適應(yīng)加密的得到網(wǎng)格,第二行是邊界條件以及自適應(yīng)網(wǎng)格下解得的速度和溫度uh,Th

3.2 實際問題

下面我們給出一個實際問題,考慮方形區(qū)域?=[0,1]×[0,1],邊界如圖2的左下角所示,其中上下邊界對溫度T和濃度C都是絕緣的,左右邊界分別給T=C=4?y?(1?y),速度在邊界上都為0.從圖2中我們可以看出在速度和溫度變化相對大的地方,網(wǎng)格能夠自適應(yīng)的進(jìn)行調(diào)整,表明改估計子對實際問題也十分有效.

參考文獻(xiàn)：

[1]Babuvˇska I,Rheinboldt W C.Error estimates for adaptive f i nite element computations[J].SIAM Journal on Numerical Analysis,1978,15(4):736-754.

[2]Eriksson K,Johnson C.Adaptive f i nite element methods for parabolic problems I:A linear model problem[J].SIAM Journal on Numerical Analysis,1991,28(1):43-77.

[3]Bangerth W,Rannacher R.Adaptive Finite Element Methods for Diあerential Equations[M].Suwitzerland:Birkh¨auser,2013.

[4]Ishak A,Nazar R,Pop I.Steady and unsteady boundary layers due to a stretching vertical sheet in a porous medium using Darcy-Brinkman equation model[J].International Journal of Applied Mechanics and Engineering,2006,11(3):623-637.

[5]Goyeau B,Songbe J P,Gobin D.Numerical study of double-diあusive natural convection in a porous cavity using the Darcy-Brinkman formulation[J].International Journal of Heat and Mass Transfer,1996,39(7):1363-1378.

[6]Liu H,Patil P R,Narusawa U.On Darcy-Brinkman equation:viscous f l ow between two parallel plates packed with regular square arrays of cylinders[J].Entropy,2007,9(3):118-131.

[7]Verf¨urth R.A Review of A Posteriori Error Estimation and Adaptive Mesh-Ref i nement Techniques[M].New York:John Wiley Sons Inc,1996.

[8]Ainsworth M,Oden J T.A Posteriori Error Estimation in Finite Element Analysis[M].New York:John Wiley Sons,2011.

[9]Verf¨urth R.A posteriori error estimates for nonlinear problems.Finite element discretizations of elliptic equations[J].Mathematics of Computation,1994,62(206):445-475.

[10]Ainsworth M,Oden J T.A posteriori error estimation in f i nite element analysis[J].Computer Methods in Applied Mechanics and Engineering,1997,142(1-2):1-88.

[11]Bochev P B,Dohrmann C R,Gunzburger M D.Stabilization of low-order mixed f i nite elements for the Stokes equations[J].SIAM Journal on Numerical Analysis,2006,44(1):82-101.

[12]Zhang Y,Hou Y,Zuo H.A posteriori error estimation and adaptive computation of conduction convection problems[J].Applied Mathematical Modelling,2011,35(5):2336-2347.