張見明+李湘賀+陸陳俊+李光耀

文章編號：16742974(2014)05005807

收稿日期：20130908

基金項目：國家自然科學(xué)基金資助項目（11172098）

作者簡介：張見明（1965-），男，湖北孝昌人，湖南大學(xué)教授，博士生導(dǎo)師

通訊聯(lián)系人，Email:zhangjm@hnu.edu.cn

摘 要：將邊界面法應(yīng)用于多域問題穩(wěn)態(tài)熱傳導(dǎo)的計算，提出了一種新的實現(xiàn)方法，稱之為多域邊界面法(MDBFM).在邊界面的基礎(chǔ)上，仿照單域問題，推導(dǎo)了多域問題穩(wěn)態(tài)熱傳導(dǎo)的矩陣組裝的方式，得出其離散的邊界積分方程，溫度等未知量即可求出.將多域邊界面法應(yīng)用于某包含62個澆筑層（共125個域）的大壩的穩(wěn)態(tài)熱傳導(dǎo)分析，得出其溫度場分布圖，并和有限元計算結(jié)果進(jìn)行了比較，溫度最大值均為20.7 ℃，且溫度分布等值線高度吻合，但多域邊界面法采用了更少的網(wǎng)格.對大壩這樣的大規(guī)模工程問題進(jìn)行計算的結(jié)果證明，本文所提出的多域邊界面法可以應(yīng)用于穩(wěn)態(tài)熱傳導(dǎo)問題，并且相較于其他方法（例如有限元法），具有同等精度，并且消耗更少的人力.

關(guān)鍵詞：邊界面法;穩(wěn)態(tài)熱傳導(dǎo);多域問題;大規(guī)模工程應(yīng)用

中圖分類號：TH164 文獻(xiàn)標(biāo)識碼：A

Application of Multidomain Boundary Face Method 

in Steady State Heat Conduction



ZHANG Jianming, LI Xianghe, LU Chenjun, LI Guangyao

(State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,

Hunan Univ, Changsha, Hunan 410082，China)

Abstract:This paper presenteda multidomain boundary face method (MDBFM) to solve steadystate heat conduction problems in largescale engineering structures. Based on BFM, and using a similar approach as that in singledomain problems, this paper derived the matrix assembly process in multidomain problems. Then, temperature and other unknown quantities can be drawn from the discrete boundary integral equation. The multidomain BFM was developed and applied to the steady state heat conduction analysis of an actual gravity dam, which contained 62 concrete layers (totally 125 domains). In order to verifythe result, the temperature contours of the dam were compared with the result obtained from FEM (using Abaqus). The comparison chart has shown that the maximum temperature value of both results is 20.7 ℃. Moreover, the contours are highly consistent. Numerical results have demonstrated that our method can achieve comparable accuracy than other methods (e.g. the FEM) at much lower cost in terms of both computer resources and human labor.

Key words:boundary face method；steady state heat conduction； multidomain problem； large scale problem



計算機(jī)輔助工程(CAE)對于推動產(chǎn)品研發(fā)具有重大的意義.目前CAE的發(fā)展已比較成熟，許多成熟的商業(yè)CAE軟件正在被廣泛地采用.但CAE仍面臨許多難題[1]，比如如何離散復(fù)雜幾何的單元才能進(jìn)行有效計算，如何處理大規(guī)模工程問題的數(shù)值計算.陸續(xù)涌現(xiàn)出的有限差分法(FDM)、有限體積法(FVM)、有限元法(FEM)和邊界元法(BEM)[2]都未能很好地解決這些問題.而目前大多數(shù)CAE軟件采用的是有限元法.

有限元法需要對整個求解域進(jìn)行離散，將會產(chǎn)生一個很大的代數(shù)方程組.對于求解三維(3D)復(fù)雜的實體，尤其是含有細(xì)小特征時，離散為可以進(jìn)行有效計算的實體單元往往比較困難.并且，有限元法的實現(xiàn)基于所求物理問題控制方程和邊界條件的等效積分“弱”形式，其試函數(shù)要求至少具有一階連續(xù)性.導(dǎo)致應(yīng)力精度總是比位移精度低一階，但在實際問題中更關(guān)注于應(yīng)力值，比如產(chǎn)生應(yīng)力集中的部位及其最大值.

相較于FEM，邊界元法(BEM)[3] 彌補(bǔ)了有限元法的不足，是一種更加有效的數(shù)值方法.它具有等幾何分析的特點，使其便于模擬復(fù)雜的幾何形狀.邊界元法基于邊界積分方程，只需對求解模型的邊界進(jìn)行離散，使求解的問題域降低了一級，大大簡化了分析和計算.并且邊界積分方程采用問題的解析基本解，具有更高的精度.然而，在傳統(tǒng)的邊界元法中，三維CAD幾何模型被離散成邊界元分析模型后，CAD模型的原始幾何信息基本被丟掉，這會從根本上導(dǎo)致計算精度的問題，有些甚至對計算起著決定性的作用，并且導(dǎo)致設(shè)計和分析成為了兩個相互獨立的過程.由于CAD模型和分析模型的分離，在邊界元自適應(yīng)網(wǎng)格細(xì)分過程中，需要反復(fù)地與CAD系統(tǒng)進(jìn)行交互，而每個交互的過程是很繁瑣的.

為了克服上述缺點，張見明教授在邊界元法[4]和邊界點法的基礎(chǔ)上，創(chuàng)造性地提出了邊界面法(BFM).在該方法中，對邊界的數(shù)值積分和場變量的插值都在邊界曲面的二維參數(shù)空間中進(jìn)行，CAE分析是直接在CAD模型上進(jìn)行的，實現(xiàn)了復(fù)雜結(jié)構(gòu)的CAE分析自動化.由于CAE模型與CAD幾何模型融為了一體，不管網(wǎng)格離散有多么粗糙，分析模型在幾何上都是精確的.并且，自適應(yīng)網(wǎng)格細(xì)分過程中,不需要再與CAD系統(tǒng)反復(fù)地進(jìn)行交互, 使自適應(yīng)分析變得簡單.

目前邊界面法已經(jīng)取得了許多研究成果，覃先云等人開發(fā)了一套參數(shù)曲面內(nèi)網(wǎng)格自動生成的算法[5]，使BFM向?qū)嶋H工程應(yīng)用邁出一大步.在覃先云的方法中，單元定義和網(wǎng)格的劃分都在邊界表面的參數(shù)空間中進(jìn)行.谷金良等人提出在BFM中用B樣條插值方案來做物理變量的近似[6-8]，并將其成功地應(yīng)用于穩(wěn)態(tài)熱傳導(dǎo)問題和彈性靜力學(xué)問題，并且通過在幾何模型中采用相同的插值方案，谷金良等人實現(xiàn)了邊界面法的等幾何分析.莊超等人實現(xiàn)了基于幾何模型的BFM[9]，其幾何模型通過細(xì)分曲面成型技術(shù)構(gòu)建.謝貴重等人通過引入改進(jìn)的距離變換、指數(shù)變換，提出了一套高效的近奇異積分計算方案[10]，解決了薄型結(jié)構(gòu)分析精度下降的問題.張見明等人通過快速多級子算法、自適應(yīng)交叉擬合（ACA）和分級矩陣（Hmatrix）技術(shù)，成功地將計算量級從O(N2)降至O(NlogN)，并且基于GPU高性能并行計算，在CUDA編程環(huán)境中實現(xiàn)邊界面法正則積分的并行加速[11]，使邊界面法的大規(guī)模工程應(yīng)用成為可能.

以上所提及的邊界面法的應(yīng)用，都是基于單域模型.然而，在工程實際分析中，往往存在許多復(fù)雜的結(jié)構(gòu)，這些結(jié)構(gòu)需要被分割成多個子域.本文將邊界面法的應(yīng)用延伸至多域模型的穩(wěn)態(tài)熱傳導(dǎo)問題.在多域問題中，矩陣的組裝是至關(guān)重要的，文中以3個兩兩相交的域為例，給出矩陣組裝的過程.然后，以某大壩為例，按照大壩的真實施工過程、施工參數(shù)和現(xiàn)場實驗所得的材料參數(shù)，并查閱資料，進(jìn)行仿真分析，研究壩體的溫度分布狀況，給出大壩采用多域邊界面法進(jìn)行穩(wěn)態(tài)熱傳導(dǎo)計算的數(shù)值結(jié)果，并和有限元法（采用Abaqus軟件）的結(jié)果進(jìn)行對比.

1 三維穩(wěn)態(tài)熱傳導(dǎo)問題的邊界面法

1.1 邊界積分方程

三維穩(wěn)態(tài)熱傳導(dǎo)問題可用如下表示：

u,ii=0,x∈Ω

u=,x∈Γu

u,ini≡q=,x∈Γq(1)

式中Ω的整個邊界為Γ.Γu和Γq分別為位勢已知邊界和法向流已知邊界，和分別是其邊界值.n是邊界外法向矢量，ni是法矢量分量，i=1,2,3. 

把上述問題轉(zhuǎn)化為等效邊界積分形式：

∫Γ(u(s)－u(y))qs(s,y)dΓ=∫Γq(s)us(s,y)dΓ (2)

式中y,s分別為邊界上的源點和場點；q=u/n代表邊界法向流量；us(s, y) 和 qs(s, y) 代表相應(yīng)的基本解，它們分別滿足：

us,ii=δ(y,s),s∈Ω(3)

qs=us/n(s) (4)

對于三維問題，

us(s,y)=14π1r(s,y) (5)

其中r是源點和場點之間的距離.

1.2 邊界積分方程的離散

本節(jié)介紹一種基于表面單元的Lagrange近似.面單元定義在表面的二維參數(shù)空間而不是在物理空間或者其他的參數(shù)空間，這樣所考慮的邊界積分的幾何數(shù)據(jù)可以通過參數(shù)轉(zhuǎn)換直接進(jìn)行計算，這樣的參數(shù)變換與參數(shù)空間的映射方案相同.換言之，可以精確得到積分中的幾何信息.下面以四節(jié)點四邊形單元為例（如圖1所示），進(jìn)行物理變量的近似.



圖1 四節(jié)點面單元及其坐標(biāo)變換

Fig.1 Fournode surface element and coordinate mapping



在這里構(gòu)建的形函數(shù)如下：

N1=14(1+ξ)(1+η), N2=14(1－ξ)(1+η)

N3=14(1－ξ)(1－η), N4=14(1+ξ)(1－η) (6)

這些形函數(shù)只用于邊界表面上的物理變量的近似，幾何數(shù)據(jù)保持精確值，這是和傳統(tǒng)邊界元法的主要區(qū)別.

u(x,y,z)=u(u,v)=u(ξ,η)=∑Nk=1Nk(ξ,η)uk

q(x,y,z)=u(u,v)=u(ξ,η)=∑Nk=1Nk(ξ,η)qk (7)

其中x=x(u, v), y=y(u, v), z=z(u, v), uk和qk分別是邊界節(jié)點上的溫度值和法向流量，N為所有插值點的個數(shù).

采用這樣的近似方案后，邊界積分方程(2)離散為：

∑Mj=1∫Γjus(s,y)∑Nk=1Nk(s)qkdΓ(s)－

∑Mj=1∫Γjqs(s,y)∑Nk=1(Nk(s)－

Nk(y))ukdΓ(s)=0(8)

把場點分布到每一個插值點后，邊界積分方程可組裝成：

Gq－Hu=0(9)

其中

Hik=∑Mj=1∫Γjqs(s,yi)(Nk(s)－Nk(yi))dΓ(s)(10)

Gik=∑Mj=1∫Γjus(s,yi)Nk(s)dΓ(s) (11)

上式中，Γj為形函數(shù)Nk(s)值不為0的邊界單元，數(shù)量記為Num.

將方程(9)進(jìn)行變換，使未知量移到左邊，已知量移到右邊，形成線性組， 

AX=b(12)

式中X是uk或qk的未知數(shù)向量.求解方程(12)，就可以得到所有節(jié)點k(k=1,2，…，N)上未知量uk或qk的值.根據(jù)節(jié)點上的值，相應(yīng)地可以求得邊界上和域內(nèi)任意一點的勢和流量值.

2 多域邊界面法在三維穩(wěn)態(tài)熱傳導(dǎo)問題中

的應(yīng)用

本章首先介紹單域問題的矩陣組裝，然后仿照單域問題和文獻(xiàn)[12]，推導(dǎo)穩(wěn)態(tài)熱傳導(dǎo)問題的多域邊界面法的邊界積分方程及其矩陣組裝.由于給定第一類（溫度、位移）或者第二類（熱流密度、面力）邊界條件時矩陣的組裝較容易，這里我們只介紹含有對流邊界的矩陣組裝方法.

2.1 單域穩(wěn)態(tài)熱問題矩陣的組裝

單域問題的邊界積分方程為：

H11H12H13H21H22H23H31H32H33u-u=G11G12G13G21G22G23G31G32G33q (13)

其中Hij,Gij分別為式（9）中的矩陣塊，上式可寫為：

H11H21H31+H12H22H32u+H13H23H33=

G11G21G31q+G12G22G32+G13G23G33 (14)

令：

=β－α (15)

其中，

β=β1β2βn,α=α10…00α2…000…αn(16)

βi=－h(huán)u

SymboleB@

，αi=－h(huán) (17)

組裝之后為：

－G11H12H13+αG13－G21H22H23+αG23－G31H32H33+αG33qu=

－H11－H21－H31+G12G22G32+G13G23G33β(18)

解方程可求出,的值通過式（15）計算得到.

2.2 多域邊界面法中矩陣的組裝

以3個兩兩相鄰的立方體為例，推導(dǎo)多域問題的邊界積分方程，如圖2所示．



圖23個兩兩相鄰立方體示意圖

Fig.2 The sketch map of 3 cubes

在這個模型中，立方體1和立方體2的邊界相交于Γ１２, 立方體1和立方體3邊界相交于Γ１3，立方體2和立方體3邊界相交于Γ23. 兩個域相交處的溫度邊界和熱流密度邊界分別為：

u

Euclid ExtrazB@

i=u

Euclid ExtrazB@

j(19)

q

Euclid ExtrazB@

i=－q

Euclid ExtrazB@

j(20)

立方體1的邊界積分方程如式(21)所示：

HddHdnHdrHd2Hd3HndHnnHnrHn2Hn3HrdHrnHrrHr2Hr3H2dH2nH2rH22H23H3dH3nH3rH32H33uu

Euclid ExtrazB@

2u

Euclid ExtrazB@

3=

GddGdnGdrGd2Gd3GndGnnGnrGn2Gn3GrdGrnGrrGr2Gr3G2dG2nG2rG22G23G3dG3nG3rG32G33qq

Euclid ExtrazB@

2q

Euclid ExtrazB@

3 (21)

其中，下標(biāo) d,n,r分別表示第一個域溫度邊界、熱流密度邊界、對流邊界所對應(yīng)的系數(shù)矩陣塊.下標(biāo)2，3表示第一個域分別與第2、第3個域相交邊界所對應(yīng)的系數(shù)矩陣塊.加入立方體1的邊界條件：

=1d，=1n，=1r (22)

并且將式(15)代入式（21），之后把未知量移到左邊，已知量移到右邊，可得：

G1dd－H1dn(－H1dr－G1drα)－H1d2－H1d3G1d2G1d3G1nd－H1nn(－H1nr－G1nrα)－H1n2－H1n3G1n2G1n3G1rd－H1rn(－H1rr－G1rrα)－H1r2－H1r3G1r2G1r3G12d－H12n(－H12r－G12rα)－H122－H123G122G123G13d－H13n(－H13r－G13rα)－H132－H133G132G133q1du1n

1r

u

Euclid ExtrazB@

12

u

Euclid ExtrazB@

13

q

Euclid ExtrazB@

12

q

Euclid ExtrazB@

13=

H1dd－G1dn－G1drH1nd－G1nn－G1nrH1rd－G1rn－G1rrH12d－G12n－G12rH13d－G13n－G13rd1n1β1=1(23)

類似立方體1，立方體2和立方體3也分別有如下重新組裝的邊界積分方程：

G211G21d－H21n(－H21r－G21rα)－H213－H211G213G2d1G2dd－H2dn(－H2dr－G2drα)－H2d3－H2d1G2d3G2n1G2nd－H2nn(－H2nr－G2nrα)－H2n3－H2n1G2n3G2r1G2rd－H2rn(－H2rr－G2rrα)－H2r3－H2r1G2r3G231G23d－H23n(－H23r－G23rα)－H233－H231G23312q2du2nr2321232=

H21d－G21n－G21rH2dd－G2dn－G2drH2nd－G2nn－G2nrH2rd－G2rn－G2rrH23d－G23n－G23rd2n2β2=2(24)

G311G312G31d－H31n(－H31r－G31rα)－H311－H312G321G322G32d－H32n(－H32r－G32rα)－H321－H322G3d1G3d2G3dd－H3dn(－H3dr－G3drα)－H3d1－H3d2G3n1G3n2G3nd－H3nn(－H3nr－G3nrα)－H3n1－H3n2G3r1G3r2G3rd－H3rn(－H3rr－G3rrα)－H3r1－H3r21323q3du3n3ru

Euclid ExtrazB@

31

u

Euclid ExtrazB@

32=

H31d－G31n－G31rH32d－G32n－G32rH3dd－G3dn－G3drH3nd－G3nn－G3nrH3rd－G3rn－G3rrd3n3β3=3(25)

再加入交界面的邊界條件，記Fkir=Hkir+Gkirα, 把式(23)~(25)整合到一起，有：

G1dd-H1dn-F1dr-H1d2-H1d3-G1d2-G1d3

G1nd-H1nn-F1nr-H1n2-H1n3-G1n2-G1n3

G1rd-H1rn-F1rr-H1r2-H1r3-G1r2-G1r3

G12d-H12n-F12r-H122-H123-G122-G123

G13d-H13n-F13r-H132-H133-G132-G133

-H211G211G21d-H21n-F21r-H213-G213

-H2d1G2d1G2dd-H2dn-F2dr-H2d3-G2d3

-H2n1G2n1G2nd-H2nn-F2nr-H2n3-G2n3

-H2r1G2r1G2rd-H2rn-F2rr-H2r3-G2r3

-H231G231G23d-H23n-F23r-H233-G213

-H311-H312G311G312G31d-H31n-F31r

-H321-H322G321G322G32d-H32n-F32r

-H3d1-H3d2G3d1G3d2G3dd-H3dn-F3dr

-H3n1-H3n2G3n1G3n2G3nd-H3nn-F3nr

-H3r1-H3r2G3r1G3r2G3rd-H3rn-F3rrq1d

u1n

1r

u

Euclid ExtrazB@

1r

u

Euclid ExtrazB@

13

q

Euclid ExtrazB@

21

q2d

u2n

2r

u

Euclid ExtrazB@

23

u

Euclid ExtrazB@

31

u

Euclid ExtrazB@

32

q3d

u3n

3r

=1

2

3(26)

求解該方程組得到邊界上節(jié)點的溫度和流量，然后利用邊界積分方程(8)，取y為任意域內(nèi)點，則得到域內(nèi)任意點的溫度.再對邊界積分方程(8)兩邊同時求y點任意方向的導(dǎo)數(shù)，即得到y點在任意方向的流量值.由于這種方向求導(dǎo)不是對形函數(shù)而是對基本解求的，不會降低擬合精度，因此使得溫度和流量的結(jié)果具有同階的計算精度.

3 數(shù)值算例

本章分析圖3所示水壩結(jié)構(gòu)在固定水溫及固定空氣溫度的邊界條件下，壩體內(nèi)溫度的分布情況.

3.1 混凝土水壩模型

如圖3所示，考慮對稱性，模型厚度取為20 m.基巖高270 m，寬471 m.基巖上的壩體高度為171 m，不均勻地分為62層（與建設(shè)過程保持一致），第一澆筑層底部與基巖接觸部分寬157 m.大壩模型左側(cè)為上游面，右側(cè)為下游面，上游面假設(shè)蓄水高度為375 m，距壩體最頂端9 m.



圖3 某重力壩段結(jié)構(gòu)示意圖

Fig.3 The crosssection of the dam



3.2 大壩的材料參數(shù)

壩段建造采用兩種不同的混凝土材料：壩體最底端兩層采用常態(tài)混凝土，基巖材料視為與第一、二澆筑層相同，壩體其他層采用碾壓混凝土.材料參數(shù)見表1．

表1 大壩材料參數(shù)

Tab.1 Properties of roller concrete and normal concrete

熱傳導(dǎo)系數(shù)/

（kJ?(m?h)-1）

比熱/

（kJ?(kg?℃)-1）

密度/

（kg?m-3）

碾壓混凝土

9.27

0.967 2

2 400

常態(tài)混凝土

8.766

0.967 2

2 450

3．3邊界條件

基巖左右兩側(cè)及底部絕熱，模型厚度方向前后面絕熱.上游面375 m高度以上、頂部表面、下游面添加空氣邊界，當(dāng)?shù)啬昶骄鶜鉁貫?0.7 ℃.上游面蓄水區(qū)添加水溫邊界條件，其中水溫邊界如下：

1）水溫在深度0~123.4 m隨深度擬合為線性變化，其擬合函數(shù)為：

Tw=20.7－0.059 157 2×h (27)

2）在距水面123.4 m以下，溫度取深水溫度13.4 ℃.

3.4 網(wǎng) 格

為做對比分析，除了使用邊界面方法分析之外，還采用基于有限元方法的商用軟件ABAQUS 11.0對此結(jié)構(gòu)進(jìn)行穩(wěn)態(tài)熱分析，由于本模型中大壩每層的厚度較薄，因此網(wǎng)格較密.在分析過程中一共使用了47 432個二次六面體單元，共計272 300個計算節(jié)點，網(wǎng)格模型如圖4所示.



圖4 有限元網(wǎng)格

Fig.4 Mesh model of FEM

在使用邊界面法進(jìn)行分析的過程中，共劃分6 494個二次單元（包括三角形單元和四邊形單元），共29 139個計算節(jié)點，網(wǎng)格模型如圖5所示.

3.5 求解結(jié)果及對比

圖6和圖7分別為ABAQUS與BFM溫度分布的計算結(jié)果對比.



圖5 BFM網(wǎng)格

Fig.5Mesh model of BFM



圖6 有限元計算結(jié)果

Fig.6Result of FEM 



圖7 BFM計算結(jié)果

Fig.7 Result of BFM

對比上面兩圖，可以看出邊界面法與有限元法分析結(jié)果基本相同，溫度最大值均為20.7 ℃，且溫度等值線高度一致.為做更加清晰詳細(xì)的對比，我們?nèi)∪鐖D8所示的兩個截面上的節(jié)點溫度，繪成溫度曲線進(jìn)行對比．



圖8截面示意圖

Fig.8 Cross sections on the model

以x軸為橫坐標(biāo)，繪出上面兩個截面上各個節(jié)點的溫度值，組成溫度曲線，如圖9所示．

x

圖9 截面1和截面2上的溫度曲線對比圖

Fig.9 The temperature curves of the two cross sections

可見，對于有限元法和邊界面法的計算結(jié)果，截面1和截面2上的溫度曲線都保持高度一致.上述兩方面的對比可說明邊界面法的計算結(jié)果與有限元基本一致.

因此，邊界面法應(yīng)用于大規(guī)模工程結(jié)構(gòu)的穩(wěn)態(tài)熱分析，具有良好的精度.而且，在邊界面法的計算中，計算時長共用了83 s，有限元法用了257 s.并且，由于邊界面法是直接基于幾何模型進(jìn)行操作，無需定義多個域的裝配以及面間的接觸等，因此邊界面法消耗較少的人力和計算機(jī)資源，具有更高的計算效率.

4 結(jié) 論 

本文用邊界面法解決了含有62個澆筑層的大壩的穩(wěn)態(tài)熱傳導(dǎo)問題.通過和有限元法的對比，證明了計算結(jié)果的正確性.說明本文所提出的多域邊界面法可以應(yīng)用在大規(guī)模工程問題中的穩(wěn)態(tài)熱分析中，并且更加便捷，具有更高的效率.

在后續(xù)工作中，我們將實現(xiàn)多域邊界面法在瞬態(tài)熱傳導(dǎo)問題中的應(yīng)用，并且考慮通過優(yōu)化矩陣的組裝方式，來進(jìn)一步提升計算效率.

參考文獻(xiàn)

［1］ ZHANG J M, TANAKA M. Fast HdBNM for largescale thermal analysis of CNTreinforced composites[J]. Computational Mechanics, 2008, 41:777-787.

［2］ 龍述堯. 邊界單元法概論[M]. 北京：中國科學(xué)文化出版社, 2002.

LONG Shuyao. The introduction of boundary element method[M].Beijing:China Science and Culture Press, 2002.(In Chinese)

［3］ HANG J M, YAO Z H, LI H. A hybrid boundary node method[J]. International Journal for Numerical Methods in Engineering, 2002, 53:751-763.

［4］ ZHANG J M, YAO Z H, TANAKA M. The meshless regular hybrid boundary node method for 2D linear elasticity[J]. Engineering Analysis with Boundary Elements, 2003, 27:259-268. 

［5］ QIN X Y, ZHANG J M, LI G Y,et al. A finite element implementation of the boundary face method for potential problems in three dimensions[J]. Engineering Analysis with Boundary Elements,2010,34:934-943.

［6］ GU J L, ZHANG J M, LI G Y. Isogeometric analysis in BIE for 3D potential problem[J]. Engineering Analysis with Boundary Elements,2012,36:858-865.

［7］ GU J L, ZHANG J M, SHENG X M, et al. Bspline approximation in boundary face method for three dimensional linear elasticity[J]. Engineering Analysis with Boundary Elements,2011,35:1159-1167.

［8］ GU J L, ZHANG J M, SHENG X M. The boundary face method with variable approximation by bspline basis function[J]. International Journal of Computational Methods,2012,9:9-18.

［9］ ZHUANG C, ZHANG J M, QIN X Y, et al. Integration of subdivision method into boundary element analysis[J]. International Journal of Computational Methods, 2012,9:19-29.

［10］XIE G Z, ZHANG J M, QIN X Y, et al. New variable transformations for evaluating nearly singular integrals in 2D boundary element method[J]. Engineering Analysis with Boundary Elements, 2011, 35:811-817.

［11］張見明，余列祥，劉路平.基于GPU加速的邊界面法正則積分的研究[J]. 湖南大學(xué)學(xué)報:自然科學(xué)版,2013,40(3):41-45.

ZHANG Jianming,YU Liexiang,LIU Luping. Research on regular integration in boundary face method based on GPUacceleration[J]. Journal of Hunan University:Natural Sciences,2013,40(3):41-45. (In Chinese)

［12］肖金生, 錢耀南, 冉一元. 柴油機(jī)活塞熱負(fù)荷的軸對稱分區(qū)邊界元分析[J]. 武漢水運工程學(xué)院學(xué)報,1986,10(2):69-78.

XIAO Jinsheng, QIAN Yaonan, RAN Yiyuan. Axisymmetric subregion boundary element analysis of thermal loading of the disel engines piston[J]. Journal of Wuhan University of Water Transportation Engineering,1986,10(2):69-78. (In Chinese)

LONG Shuyao. The introduction of boundary element method[M].Beijing:China Science and Culture Press, 2002.(In Chinese)

［3］ HANG J M, YAO Z H, LI H. A hybrid boundary node method[J]. International Journal for Numerical Methods in Engineering, 2002, 53:751-763.

［4］ ZHANG J M, YAO Z H, TANAKA M. The meshless regular hybrid boundary node method for 2D linear elasticity[J]. Engineering Analysis with Boundary Elements, 2003, 27:259-268. 

［5］ QIN X Y, ZHANG J M, LI G Y,et al. A finite element implementation of the boundary face method for potential problems in three dimensions[J]. Engineering Analysis with Boundary Elements,2010,34:934-943.

［6］ GU J L, ZHANG J M, LI G Y. Isogeometric analysis in BIE for 3D potential problem[J]. Engineering Analysis with Boundary Elements,2012,36:858-865.

［7］ GU J L, ZHANG J M, SHENG X M, et al. Bspline approximation in boundary face method for three dimensional linear elasticity[J]. Engineering Analysis with Boundary Elements,2011,35:1159-1167.

［8］ GU J L, ZHANG J M, SHENG X M. The boundary face method with variable approximation by bspline basis function[J]. International Journal of Computational Methods,2012,9:9-18.

［9］ ZHUANG C, ZHANG J M, QIN X Y, et al. Integration of subdivision method into boundary element analysis[J]. International Journal of Computational Methods, 2012,9:19-29.

［10］XIE G Z, ZHANG J M, QIN X Y, et al. New variable transformations for evaluating nearly singular integrals in 2D boundary element method[J]. Engineering Analysis with Boundary Elements, 2011, 35:811-817.

［11］張見明，余列祥，劉路平.基于GPU加速的邊界面法正則積分的研究[J]. 湖南大學(xué)學(xué)報:自然科學(xué)版,2013,40(3):41-45.

ZHANG Jianming,YU Liexiang,LIU Luping. Research on regular integration in boundary face method based on GPUacceleration[J]. Journal of Hunan University:Natural Sciences,2013,40(3):41-45. (In Chinese)

［12］肖金生, 錢耀南, 冉一元. 柴油機(jī)活塞熱負(fù)荷的軸對稱分區(qū)邊界元分析[J]. 武漢水運工程學(xué)院學(xué)報,1986,10(2):69-78.

XIAO Jinsheng, QIAN Yaonan, RAN Yiyuan. Axisymmetric subregion boundary element analysis of thermal loading of the disel engines piston[J]. Journal of Wuhan University of Water Transportation Engineering,1986,10(2):69-78. (In Chinese)

LONG Shuyao. The introduction of boundary element method[M].Beijing:China Science and Culture Press, 2002.(In Chinese)

［3］ HANG J M, YAO Z H, LI H. A hybrid boundary node method[J]. International Journal for Numerical Methods in Engineering, 2002, 53:751-763.

［4］ ZHANG J M, YAO Z H, TANAKA M. The meshless regular hybrid boundary node method for 2D linear elasticity[J]. Engineering Analysis with Boundary Elements, 2003, 27:259-268. 

［5］ QIN X Y, ZHANG J M, LI G Y,et al. A finite element implementation of the boundary face method for potential problems in three dimensions[J]. Engineering Analysis with Boundary Elements,2010,34:934-943.

［6］ GU J L, ZHANG J M, LI G Y. Isogeometric analysis in BIE for 3D potential problem[J]. Engineering Analysis with Boundary Elements,2012,36:858-865.

［7］ GU J L, ZHANG J M, SHENG X M, et al. Bspline approximation in boundary face method for three dimensional linear elasticity[J]. Engineering Analysis with Boundary Elements,2011,35:1159-1167.

［8］ GU J L, ZHANG J M, SHENG X M. The boundary face method with variable approximation by bspline basis function[J]. International Journal of Computational Methods,2012,9:9-18.

［9］ ZHUANG C, ZHANG J M, QIN X Y, et al. Integration of subdivision method into boundary element analysis[J]. International Journal of Computational Methods, 2012,9:19-29.

［10］XIE G Z, ZHANG J M, QIN X Y, et al. New variable transformations for evaluating nearly singular integrals in 2D boundary element method[J]. Engineering Analysis with Boundary Elements, 2011, 35:811-817.

［11］張見明，余列祥，劉路平.基于GPU加速的邊界面法正則積分的研究[J]. 湖南大學(xué)學(xué)報:自然科學(xué)版,2013,40(3):41-45.

ZHANG Jianming,YU Liexiang,LIU Luping. Research on regular integration in boundary face method based on GPUacceleration[J]. Journal of Hunan University:Natural Sciences,2013,40(3):41-45. (In Chinese)

［12］肖金生, 錢耀南, 冉一元. 柴油機(jī)活塞熱負(fù)荷的軸對稱分區(qū)邊界元分析[J]. 武漢水運工程學(xué)院學(xué)報,1986,10(2):69-78.

XIAO Jinsheng, QIAN Yaonan, RAN Yiyuan. Axisymmetric subregion boundary element analysis of thermal loading of the disel engines piston[J]. Journal of Wuhan University of Water Transportation Engineering,1986,10(2):69-78. (In Chinese)