多域邊界面法在穩(wěn)態(tài)熱傳導(dǎo)問題中的應(yīng)用*

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

(湖南大學(xué) 汽車車身先進設(shè)計制造國家重點實驗室,湖南 長沙 410082)

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

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

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

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

目前邊界面法已經(jīng)取得了許多研究成果，覃先云等人開發(fā)了一套參數(shù)曲面內(nèi)網(wǎng)格自動生成的算法[5]，使BFM向?qū)嶋H工程應(yīng)用邁出一大步.在覃先云的方法中，單元定義和網(wǎng)格的劃分都在邊界表面的參數(shù)空間中進行.谷金良等人提出在BFM中用B樣條插值方案來做物理變量的近似[6-8]，并將其成功地應(yīng)用于穩(wěn)態(tài)熱傳導(dǎo)問題和彈性靜力學(xué)問題，并且通過在幾何模型中采用相同的插值方案，谷金良等人實現(xiàn)了邊界面法的等幾何分析.莊超等人實現(xiàn)了基于幾何模型的BFM[9]，其幾何模型通過細分曲面成型技術(shù)構(gòu)建.謝貴重等人通過引入改進的距離變換、指數(shù)變換，提出了一套高效的近奇異積分計算方案[10]，解決了薄型結(jié)構(gòu)分析精度下降的問題.張見明等人通過快速多級子算法、自適應(yīng)交叉擬合(ACA)和分級矩陣(H-matrix)技術(shù)，成功地將計算量級從O(N2)降至O(NlogN)，并且基于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ù)，并查閱資料，進行仿真分析，研究壩體的溫度分布狀況，給出大壩采用多域邊界面法進行穩(wěn)態(tài)熱傳導(dǎo)計算的數(shù)值結(jié)果，并和有限元法(采用Abaqus軟件)的結(jié)果進行對比.

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

1.1 邊界積分方程

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

u,ii=0,?x∈Ω

(1)

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

(2)

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

us,ii=δ(y,s),?s∈Ω

(3)

qs=?us/?n(s)

(4)

對于三維問題，

(5)

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

1.2 邊界積分方程的離散

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

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

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

(6)

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

(7)

其中x=x(u,v),y=y(u,v),z=z(u,v),uk和qk分別是邊界節(jié)點上的溫度值和法向流量，N為所有插值點的個數(shù).采用這樣的近似方案后，邊界積分方程(2)離散為：

Nk(y))ukdΓ(s)=0

(8)

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

Gq-Hu=0

(9)

其中

(10)

(11)

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

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

AX=b

(12)

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

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

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

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

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

(13)

其中Hij,Gij分別為式(9)中的矩陣塊，上式可寫為：

(14)

令：

(15)

其中，

(16)

βi=-hu，αi=-h

(17)

組裝之后為：

(18)

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

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

圖2 3個兩兩相鄰立方體示意圖

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

(19)

(20)

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

(21)

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

(22)

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

(23)

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

(24)

(25)

(26)

求解該方程組得到邊界上節(jié)點的溫度和流量，然后利用邊界積分方程(8)，取y為任意域內(nèi)點，則得到域內(nèi)任意點的溫度.再對邊界積分方程(8)兩邊同時求y點任意方向的導(dǎo)數(shù)，即得到y(tǒng)點在任意方向的流量值.由于這種方向求導(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)示意圖

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

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

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

3．3 邊界條件

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

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

Tw=20.7-0.059 157 2×h

(27)

2)在距水面123.4 m以下，溫度取深水溫度13.4 ℃.

3.4 網(wǎng) 格

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

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

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

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

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

圖5 BFM網(wǎng)格

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

圖7 BFM計算結(jié)果

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

圖8 截面示意圖

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

x

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

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

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)化矩陣的組裝方式，來進一步提升計算效率.

[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,etal.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,etal.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,etal.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,etal.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]肖金生,錢耀南,冉一元.柴油機活塞熱負荷的軸對稱分區(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 engine’s piston[J].Journal of Wuhan University of Water Transportation Engineering,1986,10(2):69-78.(In Chinese)