虛邊界元法在二維涂層結(jié)構(gòu)溫度場(chǎng)中的應(yīng)用

王發(fā)杰， 張耀明

(山東理工大學(xué) 理學(xué)院， 山東 淄博 255091)

近年來(lái)，隨著表面技術(shù)及工程的發(fā)展，各種功能涂層由于其具有眾多良好的性能而愈來(lái)愈引起人們的重視，應(yīng)用范圍涉及能源、石油化工、建筑、機(jī)械、航空航天等諸多領(lǐng)域，與國(guó)家經(jīng)濟(jì)建設(shè)、國(guó)防及人們的日常生活關(guān)系日益緊密，已成為表面工程技術(shù)的一個(gè)重要組成部分[1].然而，涂層材料厚度一般較薄，約在微米級(jí)甚至納米級(jí)，受其厚度尺寸的限制，涂層材料中物理量的數(shù)值分析一直是工程中的難點(diǎn).對(duì)這類結(jié)構(gòu)采用有限元分析，受結(jié)構(gòu)形狀制約，為保證單元不退化，其長(zhǎng)寬比必須協(xié)調(diào)，故所需劃分的單元數(shù)量過(guò)大[2].邊界元法可有效地處理涂層結(jié)構(gòu)問(wèn)題[3-5]，但涉及奇異邊界積分和幾乎奇異邊界積分的計(jì)算麻煩，邊界元法的精確度取決于這些奇異積分的計(jì)算效果，因此具有一定的難度而且耗時(shí).

與傳統(tǒng)邊界元法相比，虛邊界元法(Virtual boundary element method, VBEM)是一種無(wú)奇異邊界元法[6]，其基本原理是通過(guò)真實(shí)邊界外虛擬分布的源對(duì)場(chǎng)點(diǎn)的影響疊加產(chǎn)生一個(gè)位勢(shì)積分，然后利用實(shí)邊界條件，建立虛邊界積分方程.虛邊界元法避免了奇異積分的麻煩處理，消除了傳統(tǒng)邊界元法的邊界層效應(yīng)，具有程序設(shè)計(jì)簡(jiǎn)單，精度高，耗時(shí)短等優(yōu)點(diǎn)，已廣泛應(yīng)用于各種常規(guī)結(jié)構(gòu)，并取得了很好的效果[7-9].但是用于涂層結(jié)構(gòu)的研究至今尚未發(fā)現(xiàn).本文將虛邊界元法應(yīng)用于二維涂層結(jié)構(gòu)溫度場(chǎng)問(wèn)題.由于虛邊界分布在物理區(qū)域外，所以對(duì)于涂層問(wèn)題的虛邊界元法來(lái)說(shuō)，一個(gè)區(qū)域的虛邊界可能劃分在另外一個(gè)區(qū)域.對(duì)此，本文利用復(fù)雜的多域技術(shù)(MDT)，發(fā)展了多域虛邊界元法(MD-VBEM)，有效地計(jì)算了涂層問(wèn)題.從而為涂層結(jié)構(gòu)的研究開(kāi)辟了新的途徑，同時(shí)也拓展了虛邊界元法的應(yīng)用領(lǐng)域，即使涂層結(jié)構(gòu)的狹長(zhǎng)比小到1E-10，依然可獲得非常高精度的數(shù)值解，而且方法程序設(shè)計(jì)簡(jiǎn)單，效率較高.

1 二維薄體結(jié)構(gòu)位勢(shì)問(wèn)題的虛邊界元法

假定 Ω是R2中的一個(gè)有界區(qū)域，Γ=?Ω是邊界.n=(n1,n2)是區(qū)域Ω的邊界Γ在x點(diǎn)處的單位外法向量.

1.1 位勢(shì)問(wèn)題的控制微分方程

二維位勢(shì)問(wèn)題的控制微分方程為

邊界條件為

式中u為勢(shì)函數(shù)；n為邊界外法線；Γu、Γq分別是u和?u/?n已知的邊界.

二維位勢(shì)問(wèn)題控制方程的基本解為

1.2 二維位勢(shì)問(wèn)題虛邊界元法的積分方程

設(shè)想假定Ω被嵌入一個(gè)無(wú)限的區(qū)域中，在Γ外的無(wú)限區(qū)域中有一延拓邊界Γ′(這里稱為虛邊界)，沿著邊界Γ′有一個(gè)待定的虛擬密度函數(shù)Φ(y)，令此虛擬密度函數(shù)對(duì)真實(shí)邊界所產(chǎn)生的位勢(shì)或法向梯度與邊界條件一致，從而達(dá)到求解待定密度的目的.以上稱之為虛邊界元法.

位勢(shì)問(wèn)題中的虛邊界積分方程為

(1)

(2)

計(jì)算內(nèi)點(diǎn)位勢(shì)和位勢(shì)梯度的邊界積分方程表示為

(3)

i=1,2,x∈Ω

(4)

1.3 二維涂層結(jié)構(gòu)溫度場(chǎng)問(wèn)題的虛邊界元法

圖1 分域法結(jié)構(gòu)圖

對(duì)于虛邊界元法，在Ω1上可建立如下矩陣方程

(5)

同理，在Ω2上可建立如下矩陣方程

(6)

對(duì)于適定的邊值問(wèn)題，或者邊界上的溫度已知或者溫度梯度已知.邊界離散化后，每個(gè)節(jié)點(diǎn)上都會(huì)產(chǎn)生一個(gè)代數(shù)方程，方程的個(gè)數(shù)與虛邊界節(jié)點(diǎn)處待求密度函數(shù)的個(gè)數(shù)相同，因而可以數(shù)值求解.分域法將區(qū)域Ω1與Ω2看成兩個(gè)獨(dú)立的問(wèn)題來(lái)處理，在各自區(qū)域上利用虛邊界元法進(jìn)行計(jì)算，但在Ω1與Ω2的共同邊界ΓI上，溫度與溫度梯度都是未知的，因此未知參量的個(gè)數(shù)此時(shí)大于代數(shù)方程的個(gè)數(shù).要使得邊值問(wèn)題可解，必須引入邊界ΓI上的溫度和熱通量協(xié)調(diào)條件：

(7)

假設(shè)邊界Γ1，Γ2上節(jié)點(diǎn)的位勢(shì)已知，根據(jù)條件式(7)、式(5)和式(6)可合并成

(8)

式(8)即為涂層結(jié)構(gòu)溫度場(chǎng)虛邊界元法的基本列式.通過(guò)式(8)，可求出Ω1與Ω2虛邊界上的節(jié)點(diǎn)密度函數(shù)，進(jìn)而可以利用內(nèi)點(diǎn)積分方程求出內(nèi)點(diǎn)的物理參量.類似地，可寫出其它邊界條件下相應(yīng)的方程組.

顯然，以上過(guò)程可以直接推廣到多涂層結(jié)構(gòu)問(wèn)題，只是聯(lián)立方程的個(gè)數(shù)有所增加，這里就不再過(guò)多闡述.

2 數(shù)值算例

利用兩個(gè)涂層問(wèn)題的數(shù)值算例，來(lái)驗(yàn)證本文方法的有效性，所有算例均采用常單元等額配點(diǎn)法.為了表明方法數(shù)值解的準(zhǔn)確性，定義如下平均相對(duì)誤差

(9)

虛邊界元法的求解精度受虛實(shí)邊界間的距離影響.以下算例中，選擇近、中、遠(yuǎn)三種不同的虛實(shí)邊界間的距離[10]，對(duì)不同狹長(zhǎng)比的涂層問(wèn)題進(jìn)行計(jì)算，都獲得了令人滿意的數(shù)值結(jié)果，表明虛實(shí)邊界間的有效距離選取范圍非常寬泛.

例1 圓環(huán)涂層結(jié)構(gòu)的熱流問(wèn)題.基體Ω1內(nèi)徑為r1=10m，外徑為r2=11m；涂層Ω2外徑為r3，如圖2(a)所示.已知基體內(nèi)表面溫度為10℃，涂層外表面溫度為20℃.基體導(dǎo)熱率為k1=1，涂層的導(dǎo)熱率為k2=2.定義薄體區(qū)域特征值最小尺寸與最大尺寸之比δ=(r3-r2)/r1為狹長(zhǎng)比.

根據(jù)熱傳導(dǎo)理論，涂層與基體溫度的解析解分別為

(a) (b) (c)圖2 涂層結(jié)構(gòu)圓環(huán)區(qū)域的熱流

單元?jiǎng)澐智闆r不變，狹長(zhǎng)比為δ=10-10，虛實(shí)邊界距離取中等距離組d1=5,d2=10.圖5和圖6分別列出了不同狹長(zhǎng)比下，界面點(diǎn) B(r2,0)上溫度和熱流的計(jì)算結(jié)果.可以看出數(shù)值解與精確解十分地接近.此算例充分體現(xiàn)了本文方法在計(jì)算涂層結(jié)構(gòu)問(wèn)題時(shí)的可靠性和精確度.

圖3 界面溫度解的平均相對(duì)誤差

圖4 涂層外表面熱流解的平均相對(duì)誤差

圖5 界面點(diǎn)B處的溫度

圖6 界面點(diǎn)B處的熱流

例2 矩形涂層結(jié)構(gòu)的熱流問(wèn)題.涂層的厚度為h，基體的厚度為H=1m，寬度為 L=2m，涂層與基體的導(dǎo)熱率分別為 k1=1，k2=2，邊界條件如圖7(a)所示.

(a) (b)圖7 涂層結(jié)構(gòu)矩形區(qū)域的熱流

采用MD-VBEM，圖7(b)給出了基體和涂層結(jié)構(gòu)及虛邊界計(jì)算模型.在這里將基體和涂層的虛、實(shí)邊界距離取為相同的值d.基體虛邊界四邊各劃分為10個(gè)單元，涂層上下虛邊界各劃分為10個(gè)單元，左右兩邊各劃分為2個(gè)單元，總共64個(gè)單元.當(dāng)涂層厚度h從10-1m到10-10m之間變化時(shí)，分別取近中遠(yuǎn)三種不同的虛實(shí)邊界距離，計(jì)算界面上配點(diǎn)處的溫度以及涂層上表面配點(diǎn)處的熱流.圖8及圖9給出了數(shù)值解的平均相對(duì)誤差隨涂層厚度變化的曲線.從圖8和圖9中可以看出，當(dāng)涂層厚度h從10-1m到10-10m變化時(shí)，中等距離d=10和較遠(yuǎn)距離組d=20時(shí)的數(shù)值解相對(duì)誤差非常小；較小距離d=0.5時(shí)的數(shù)值解的精度稍差，但仍然具有較高精度.表明，雖然虛實(shí)邊界間的距離選取對(duì)解的精度有影響，但是有效距離的選取范圍任然相當(dāng)?shù)貙挿?

圖8 界面溫度解的平均相對(duì)誤差

圖9 涂層上表面熱流解的平均相對(duì)誤差

單元?jiǎng)澐智闆r不變，涂層厚度為h=1.0E-10m，虛實(shí)邊界距離取中等距離d=10.圖10、圖11分別給出了界面點(diǎn)B(0.5,1)處熱流和溫度梯度 ?u/?x1的計(jì)算結(jié)果.從圖10、圖11中可以看出，所得數(shù)值解與精確解十分吻合，并不受涂層厚度變化的影響，即使涂層厚度h達(dá)到1.0E-10m，本文方法任然可以得到可靠、穩(wěn)定的數(shù)值解.

圖10 界面點(diǎn)B(0.5,1)處的熱流

圖11 界面點(diǎn) B(0.5,1)處的溫度梯度

3 結(jié)束語(yǔ)

本文將虛邊界元法應(yīng)用于二維涂層結(jié)構(gòu)溫度場(chǎng)問(wèn)題，提出了多域虛邊界元法，成功求解了二維涂層結(jié)構(gòu)問(wèn)題.從而給出了求解二維涂層結(jié)構(gòu)溫度場(chǎng)問(wèn)題的新途徑，同時(shí)也拓展了虛邊界元法的應(yīng)用范圍.數(shù)值算例表明，虛實(shí)邊界距離的選取相當(dāng)寬泛，即使涂層結(jié)構(gòu)的狹長(zhǎng)比小到10-10，依然可獲得高精度的數(shù)值解.

[1] 胡傳炘. 特種功能涂層[M].北京: 北京工業(yè)大學(xué)出版社, 2009.

[2] Luo J F, Liu Y J, Berger E J. Analysis of two-dimensional thin structures (from micro- to nano-scales) using the boundary element method[J].Computational Mechanics, 1998, 22:404-412.

[3] Du F, Lovell M R, Wu T W. Boundary element method analysis of temperature fields in coated cutting tools[J]．International Journal of Solids and Structures, 2001, 38:4 557-4 570.

[4] 程長(zhǎng)征,牛忠榮,周煥林,等.涂層結(jié)構(gòu)中溫度場(chǎng)的邊界元分析[J].合肥工業(yè)大學(xué)學(xué)報(bào), 2006, 29(3): 326-329.

[5] 張耀明,谷巖.涂層結(jié)構(gòu)中溫度場(chǎng)的邊界元解[J].固體力學(xué)學(xué)報(bào), 2011，32(2)：133-141.

[6] 孫煥純，許強(qiáng). 無(wú)奇異邊界元法[M].大連:大連理工大學(xué)出版社, 1999.

[7] Yao W A, Wang H, Virtual boundary element method for 2-D piezpelectric media[J]. Finite Elements Anal Des,2005,41:875-891.

[8] Li X C, Yao W A.Virtual boundary element-integral collocation method for the plane magnetoelectroelastic solids Engineering[J]. Analysis with Boundary Elements,2006,30:709-717.

[9] Yang D S, Xu Q, Virtual boundary meshless least square integral method with moving least squares approximation for 2D elastic problem[J]. Engineering Analysis with Boundary Elements 2013,37:616-623.

[10] 張耀明,孫煥純,楊家新.虛邊界元法的理論分析[J].計(jì)算力學(xué)學(xué)報(bào), 2000, 17(1): 56-62.