王發(fā)杰， 張耀明

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

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

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

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

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

1.1 位勢問題的控制微分方程

二維位勢問題的控制微分方程為

邊界條件為

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

二維位勢問題控制方程的基本解為

1.2 二維位勢問題虛邊界元法的積分方程

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

位勢問題中的虛邊界積分方程為

(1)

(2)

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

(3)

i=1,2,x∈Ω

(4)

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

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

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

(5)

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

(6)

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

(7)

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

(8)

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

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

2 數(shù)值算例

利用兩個涂層問題的數(shù)值算例，來驗證本文方法的有效性，所有算例均采用常單元等額配點法.為了表明方法數(shù)值解的準確性，定義如下平均相對誤差

(9)

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

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

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

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

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

圖3 界面溫度解的平均相對誤差

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

圖5 界面點B處的溫度

圖6 界面點B處的熱流

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

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

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

圖8 界面溫度解的平均相對誤差

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

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

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

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

3 結(jié)束語

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

[1] 胡傳炘. 特種功能涂層[M].北京: 北京工業(yè)大學出版社, 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] 程長征,牛忠榮,周煥林,等.涂層結(jié)構(gòu)中溫度場的邊界元分析[J].合肥工業(yè)大學學報, 2006, 29(3): 326-329.

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

[6] 孫煥純，許強. 無奇異邊界元法[M].大連:大連理工大學出版社, 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].計算力學學報, 2000, 17(1): 56-62.