孫蘭銀,蘇芳明

(信陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 河南 信陽 464000)

0 引言

在自然界中,熱傳導(dǎo)、電磁勢和流體等問題都可用偏微分方程來描述,它們的解一般很難用解析公式表達(dá)出來,運(yùn)用數(shù)值方法來計(jì)算它們的近似解是十分重要的。有限元方法是現(xiàn)代科學(xué)與工程計(jì)算領(lǐng)域中求解偏微分方程的一個(gè)主要工具。1943年,數(shù)學(xué)家COURANT[1]提出了有限元方法；20世紀(jì)50年代,有限元方法是作為解決固體力學(xué)問題的方法出現(xiàn)的；1960年,CLOUGH[2]在處理平面彈性問題時(shí),第一次提出了“有限元方法”這一名稱,并成功應(yīng)用于飛機(jī)結(jié)構(gòu)的分析；20世紀(jì)60年代,馮康發(fā)表的論文《基于變分原理的差分格式》[3],標(biāo)志著有限元方法在我國的誕生?；趥鹘y(tǒng)的有限元方法,學(xué)者們也提出了許多新的數(shù)值計(jì)算方法,例如:混合有限元方法[4]、自適應(yīng)有限元法[5]、DG方法[6]、弱Galerkin有限元法[7]。

有限元基函數(shù)的選取會(huì)對數(shù)值解的精確度產(chǎn)生影響。為解決高階多項(xiàng)式曲線的表示問題,文獻(xiàn)[8-9]等對混合代數(shù)和三角多項(xiàng)式空間進(jìn)行了擴(kuò)展,在空間T=span{1,t,…,tp-2,sint,cost}上定義了p次C-Bézier基函數(shù)以及空間T=span{sinht,cosht,tq-2,…,t,1}上的q次H-Bézier基函數(shù)。

如何得到一個(gè)精確度高且計(jì)算量較小的方法一直是有限元方法研究中的一個(gè)重要課題。HUGES等[10]將NURBS基函數(shù)應(yīng)用于有限元分析，BHATTI等[11]提出在求解偏微分方程中使用Bernstein基函數(shù)。本文考慮一維熱傳導(dǎo)方程初邊值問題的數(shù)學(xué)模型,在有限元方法中利用C-Bézier和H-Bézier基函數(shù)構(gòu)造試探函數(shù)和測試函數(shù)空間,通過數(shù)值實(shí)驗(yàn)來驗(yàn)證理論結(jié)果的正確性與方法的可行性。

1 一維熱傳導(dǎo)方程

熱傳導(dǎo)問題是描述在某個(gè)區(qū)域內(nèi)物體的溫度分布規(guī)律。給定一根均勻且各向同性的細(xì)桿,長為L(L≥0),假設(shè)細(xì)桿內(nèi)部有熱源且與周邊介質(zhì)有熱交換,f0(x,t)為熱源強(qiáng)度,C為熱容,ρ為密度,k為傳導(dǎo)率,u=u(x,t)描述在t時(shí)刻x位置材料的溫度,則Dirichlet邊界條件下的一維熱傳導(dǎo)方程初邊值問題為[12-13]：

(1)

一維熱傳導(dǎo)方程初邊值問題的變分形式為:求u=u(x,t)∈H1(Ω)(t為參數(shù)),滿足

(2)

在H1(Ω)中取一n維子空間Vh,令Vh0表示Vh中在邊界處值為0的函數(shù)構(gòu)成的空間,則半離散的Gakerkin方法為：求uh(x,t)∈Vh,滿足

(3)

進(jìn)一步利用差分法對時(shí)間t離散,得到一維熱傳導(dǎo)方程的全離散化Galerkin方法。

2 C-Bézier和H-Bézier基函數(shù)

(4)

其中:p≥2,ω∈[0,α],α∈(0,π],

α為形狀參數(shù)[8]。

(5)

其中:q≥2,ω∈[0,α],α∈(0,+∞),

α為形狀參數(shù)[9]。

C-Bézier和H-Bézier基函數(shù)都具有單位分解、端點(diǎn)插值、升階等性質(zhì),具體性質(zhì)及證明見參考文獻(xiàn)[8-9]。

定義3 稱參數(shù)曲線段

(6)

圖1 C-Bézier曲線Fig.1 C-Bézier curves

定義4 稱參數(shù)曲線段

(7)

圖2 H-Bézier曲線Fig.2 H-Bézier curves

C-Bézier和H-Bézier基函數(shù)引入了形狀參數(shù)α,α可以調(diào)節(jié)曲線的形狀和幾何性質(zhì)。以C-Bézier基函數(shù)為例,調(diào)節(jié)區(qū)域?yàn)锽ézier曲線到α=π時(shí)的C-Bézier曲線之間。曲線首端和末端的位置及切矢方向不隨形狀參數(shù)α的改變而改變。同時(shí)在多段C-Bézier曲線進(jìn)行G1拼接時(shí),α改變不會(huì)影響曲線之間原有的連續(xù)性,為曲線曲面造型帶來了便利性和靈活性。

3 以C-Bézier和H-Bézier為基函數(shù)的有限元方法

為C-Bézier或H-Bézier基函數(shù)構(gòu)成的測試函數(shù)空間。

當(dāng)p=2時(shí),二次C-Bézier基函數(shù)為:

(8)

其中ω∈[0,α],α∈(0,π]。

當(dāng)q=2時(shí),二次H-Bézier基函數(shù)為:

(9)

圖3 二次C-Bézier和H-Bézier基函數(shù)Fig. 3 Quadratic C-Bézier and H-Bézier basis

(10)

其中x∈[xn,xn+1],α∈(0,π]。

(11)

用向量表示為：

MX′(t)+A(t)X(t)=b(t),

(12)

進(jìn)一步用差分法對時(shí)間t離散,得到全離散的Galerkin方法。

4 穩(wěn)定性及誤差分析

全離散格式是由常微分方程組(11)離散化得到的,因此需要討論它的穩(wěn)定性。使用Crank-Nicolson格式對時(shí)間t進(jìn)行離散:

(f,vh),

其中:上標(biāo)n表示在t=tn=nτ處的近似。取

得

利用a(·,·)的對稱性得

由Poincaré不等式可以證得穩(wěn)定性。

由Céa引理[14]可得,有限元的誤差估計(jì)轉(zhuǎn)化為插值誤差估計(jì)。當(dāng)一維熱傳導(dǎo)方程半離散時(shí),即僅對空間離散,L2范數(shù)誤差估計(jì)為

‖u-uh‖≤C1h2‖u‖2,

其中：C1為常數(shù),h為步長,利用Aubin-Nitche技巧可證明之。

引理1 設(shè)u和uh分別為方程(1)和(3)的解,對任意的t∈[0,+∞),誤差估計(jì)為

‖uh(t)-u(t)‖≤

其中：C2為常數(shù)。

引理1的證明參見文獻(xiàn)[14],對時(shí)間t進(jìn)行離散,由引理1可以得到一維熱傳導(dǎo)方程全離散化的誤差估計(jì)。

5 數(shù)值算例

運(yùn)用Lagrange、C-Bézier和H-Bézier基函數(shù)的有限元法求解熱傳導(dǎo)方程并比較數(shù)值解的精度。在下述算例中,均使用二次的基函數(shù)構(gòu)成有限元方法中的試探函數(shù)與測試函數(shù)空間,有限元節(jié)點(diǎn)相同,所用網(wǎng)格一致。

例1 考慮一根均勻且長度為1 m的細(xì)桿,初始溫度分布為0 ℃,熱擴(kuò)散系數(shù)c=2,熱容為1,細(xì)桿內(nèi)部有熱源,在t時(shí)刻,熱源強(qiáng)度f0(x,t)=xcos(xt)+2t2sin(xt),細(xì)桿左端保持溫度0 ℃,右端溫度為sint℃,運(yùn)用有限元方法求解細(xì)桿內(nèi)部溫度分布,并分析誤差。

其一維熱傳導(dǎo)方程的初邊值問題為：

(13)

表1 方程(13)在t=2時(shí)刻節(jié)點(diǎn)處的無窮范誤差Tab. 1 Infinite norm errors at nodes of Eq.(13) for t=2

圖4 當(dāng) 時(shí), 方程(13)誤差圖像(a)和細(xì)桿內(nèi)部溫度隨時(shí)間分布圖(b)Fig.4 The error graph (a) at t=2 and interior temperature distribution diagram (b) of thin rod of Eq. (13) over time

其一維熱傳導(dǎo)方程邊值問題為

(14)

表2 方程(14)在t=1時(shí)刻在節(jié)點(diǎn)處的無窮范誤差Tab. 2 Infinite norm errors at nodes of Eq.(14) for t=1

圖5 當(dāng)時(shí), 方程(14)誤差圖像(a)和細(xì)桿內(nèi)部溫度隨時(shí)間分布圖(b)Fig. 5 The error graph (a) at t=1 and interior temperature distribution diagram (b) of thin rod of Eq. (14) over time

例3 考慮一根均勻且長為1 m的細(xì)桿, 初始溫度分布為0 ℃, 熱擴(kuò)散系數(shù)c和熱容都為1。 細(xì)桿內(nèi)部有熱源,在t時(shí)刻, 熱源強(qiáng)度f0(x,t)=tsin(πx)(2+π2t), 細(xì)桿兩端溫度為0 ℃。 運(yùn)用有限元方法求解細(xì)桿內(nèi)部溫度分布, 并分析誤差。 其一維熱傳導(dǎo)方程的初邊值問題為：

(15)

表3 方程(15)在t=1時(shí)刻在節(jié)點(diǎn)處的無窮范誤差Tab. 3 Infinite norm errors at nodes of Eq.(15) for t=1

圖6 當(dāng)時(shí), 方程(15)誤差圖像(a)和細(xì)桿內(nèi)部溫度隨時(shí)間分布圖(b)Fig. 6 The error graph (a) at t=1 and interior temperature distribution diagram (b) of thin rod of Eq. (15) over time

由例3可見, 運(yùn)用C-Bézier基函數(shù)構(gòu)成測試函數(shù)及試探函數(shù)空間, 在步長不變的情況下, 通過選取適當(dāng)?shù)男螤顓?shù)α, 此時(shí)α=πh, 得到的數(shù)值解的精度比Lagrange型基函數(shù)高三個(gè)數(shù)量級以上, 說明C-Bézier基函數(shù)的有限元法對熱傳導(dǎo)過程的模擬效果更好。

6 總結(jié)

介紹了C-Bézier和H-Bézier基函數(shù)的定義及良好性質(zhì), 并用其構(gòu)造有限元方法中的測試函數(shù)與試探函數(shù)空間, 通過例子說明利用這兩種基函數(shù)求出的數(shù)值解與精確解的誤差更小, 精度更高, 證明了理論結(jié)果。

形狀參數(shù)α的選取會(huì)影響數(shù)值解的精確度, 同時(shí)在進(jìn)行誤差估計(jì)時(shí), 常數(shù)C1、C2的確定也非常重要。 下一步將繼續(xù)研究如何選取最優(yōu)的形狀參數(shù)α, 使得數(shù)值解與精確解之間的誤差達(dá)到最小, 并探討如何利用簡單的方法確定常數(shù)C1、C2進(jìn)行誤差估計(jì)。