高 宇， 劉 彪， 李通盛， 楊新貴，何錢生， 王 橋， 周 偉， 曹 悅

(1. 大唐宣威水電開發(fā)有限公司，云南 宣威 655400；2. 武漢大學(xué) 水利水電學(xué)院，湖北 武漢 430072)

熱傳導(dǎo)，是介質(zhì)內(nèi)無宏觀運(yùn)動(dòng)時(shí)的傳熱現(xiàn)象，無論在生活中或工程中，都十分常見，簡而言之，只要介質(zhì)內(nèi)或者介質(zhì)之間存在溫度差，就一定會(huì)發(fā)生傳熱。此外，在工程結(jié)構(gòu)中，由于各部分不能自由伸縮的結(jié)構(gòu)出現(xiàn)溫度變化或者內(nèi)部各部分溫度不同而產(chǎn)生的溫度應(yīng)力是結(jié)構(gòu)安全性的關(guān)鍵性影響因素。因此研究結(jié)構(gòu)內(nèi)部的溫度分布是十分必要的。進(jìn)一步的，混凝土結(jié)構(gòu)在澆筑初期會(huì)由于水化熱而產(chǎn)生大量熱量，或者內(nèi)部發(fā)生火災(zāi)釋放大量熱的建筑結(jié)構(gòu)，可以概括為含熱源的工程結(jié)構(gòu)，由于熱源的存在，結(jié)構(gòu)內(nèi)部溫度場(chǎng)的分布會(huì)發(fā)生較大變化，導(dǎo)致溫度應(yīng)力急劇增大或者分布發(fā)生較大變化，進(jìn)而對(duì)結(jié)構(gòu)安全性帶來巨大的安全隱患。

目前研究熱傳導(dǎo)問題的數(shù)值方法很多，例如無網(wǎng)格法[4]，有限元法[5]，邊界點(diǎn)法[6]等。邊界元法[7,8]作為一種重要的半解析數(shù)值方法，可以對(duì)研究問題進(jìn)行降維，只需要對(duì)研究域的邊界進(jìn)行離散，無需在內(nèi)部劃分網(wǎng)格，這是邊界元法十分顯著的優(yōu)勢(shì)。尤其是在研究大尺寸，無限域問題時(shí)，其優(yōu)勢(shì)更加明顯。目前，邊界元法在位勢(shì)問題、聲學(xué)、波的傳播、靜力學(xué)、動(dòng)力學(xué)等方面都已經(jīng)得到了廣泛應(yīng)用。對(duì)于一般熱傳導(dǎo)問題，邊界元法也能非常精確高效地解決。

而當(dāng)面臨考慮熱源的熱傳導(dǎo)問題時(shí)，在邊界積分方程中，會(huì)出現(xiàn)域積分，如果不加以處理，則需要進(jìn)行單元離散，使得邊界元喪失了降維的優(yōu)勢(shì)。因此，尋找一種能夠?qū)⒂蚍e分轉(zhuǎn)化為邊界積分的方法是十分必要的。

現(xiàn)存的域積分轉(zhuǎn)化方法中，應(yīng)用最為廣泛的方法是雙互易方法(Dual Reciprocity Method，DRM)[11,12]。DRM不需要進(jìn)行研究域內(nèi)的單元離散化，而是采用內(nèi)部節(jié)點(diǎn)和一系列指定的徑向基函數(shù)(Radial Basis Function，RBF)進(jìn)行逼近。DRM的精度和效率在很大程度上取決于RBF的分布和位置，內(nèi)部點(diǎn)的位置可以任意選擇，但域內(nèi)用于插值的形函數(shù)應(yīng)具有特殊解，不能任意選擇，對(duì)于復(fù)雜的DRM問題，要獲得特殊解不是一件容易的事情。此外，還有用于求解泊松方程的蒙特卡羅積分法(Monte Carlo Integration Method，MCM)[13]。此方法簡單易行，但不如DRM精確。近年來，還出現(xiàn)了一種十分強(qiáng)大的方法，稱為徑向積分法(Radial Integration Method，RIM)[14]，其可以將域積分轉(zhuǎn)化為邊界積分和徑向積分。一般情況下，徑向積分可采用解析或數(shù)值計(jì)算。二維和三維域積分的計(jì)算方法是統(tǒng)一的，也不需要對(duì)區(qū)域進(jìn)行離散。

本文采用了直線積分法(Line Integration Method，LIM)[15]處理域積分。直線積分法基于散度定理將域積分轉(zhuǎn)化為包含一維線積分的邊界積分，無需對(duì)研究域內(nèi)進(jìn)行離散。目前已經(jīng)應(yīng)用于彈性力學(xué)和位勢(shì)問題中的域積分處理，其有效性和精確度也得到了驗(yàn)證。此外，該方法還十分容易和快速多級(jí)算法進(jìn)行結(jié)合以適用于計(jì)算大規(guī)模工程問題。

1 考慮熱源的邊界積分方程

結(jié)合前人研究，可知控制方程為：

(1)

式中：k為熱傳導(dǎo)系數(shù)；B(x)為熱源項(xiàng)；θ為溫度；x為研究域Ω內(nèi)的點(diǎn)，其組成元素為xi。

邊界情況可設(shè)為：

(2)

通過引入權(quán)函數(shù)、分部積分法和散度定理，控制方程可寫成下面的積分方程：

(3)

(4)

(5)

式中：n為邊界Γ的外法向量，其組成元素為ni；r為點(diǎn)x和y的距離，即

r=‖x-y‖

(6)

可以看到，式中含有一項(xiàng)域積分：

(7)

需要利用直線積分法將域積分轉(zhuǎn)化為邊界積分。

2 直線積分法

本節(jié)將介紹利用直線積分法將二維問題中的域積分轉(zhuǎn)化為邊界積分，為了將問題一般化，現(xiàn)假設(shè)一個(gè)域積分為：

(8)

假設(shè)函數(shù)f(x,y)可以從矢量場(chǎng)的散度得到，具體表達(dá)為：

f(x,y)=divF=?·F

(9)

式中：div為散度運(yùn)算符；F為向量場(chǎng)，可表達(dá)為：

F=f1(x,y)e1+f2(x,y)e2

(10)

式中：ei為直角坐標(biāo)系中的一組基向量。

通過運(yùn)用散度定理，可以將域積分轉(zhuǎn)化為以下形式：

(11)

為了尋找一個(gè)合適的函數(shù)F滿足公式，可以采用特解法，故而式可以改寫為：

(12)

或

(13)

式中：nx，ny為向量n在x，y軸方向的分量，并且有

(14)

(15)

通過使用式，可以將域積分式轉(zhuǎn)化為包含直線積分的邊界積分，為了簡單起見，可以將任意積分線的起始點(diǎn)設(shè)為一個(gè)常數(shù)，x=a，其中a為任意值。并且可以將函數(shù)f(x,y)的定義擴(kuò)展至R2，這樣就不需要去對(duì)積分點(diǎn)在積分線上的位置做判斷。

假設(shè)將邊界劃分為N個(gè)邊界單元，這樣式可改寫為：

(16)

(17)

式(17)的計(jì)算一般都可以采用數(shù)值計(jì)算的方式，現(xiàn)在可以將計(jì)算域用背景網(wǎng)格覆蓋并逐級(jí)等分為小網(wǎng)格，進(jìn)而可以將積分線截?cái)喑蔀橐欢味蔚姆e分線，則式(17)可以改寫為：

(18)

式中：Lk為背景網(wǎng)格劃分后的直線積分線段。

引入背景網(wǎng)格之后，利用每個(gè)網(wǎng)格內(nèi)的積分線段上的高斯積分點(diǎn)就可以直接進(jìn)行計(jì)算并且可以獲得很高的精度。進(jìn)一步地，域積分可轉(zhuǎn)化為：

(19)

最后，將直線積分法應(yīng)用于邊界積分方程中的域積分的計(jì)算，具體轉(zhuǎn)換后的表達(dá)式為：

(20)

3 數(shù)值驗(yàn)證

為了驗(yàn)證直線積分法在邊界元法中域積分轉(zhuǎn)換的有效性和精確度，本文采用了兩個(gè)例子進(jìn)行驗(yàn)證，第一個(gè)例子為一個(gè)承受溫度應(yīng)力的矩形梁，并已知其解析解，進(jìn)而可以驗(yàn)證本文所采用方法的精度。第二個(gè)例子為一個(gè)混凝土重力壩，計(jì)算結(jié)果將與有限元法模擬結(jié)果進(jìn)行對(duì)比以驗(yàn)證本方法在大壩溫度場(chǎng)分析中的有效性。

3.1 矩形柱的含熱源熱傳導(dǎo)分析

如圖1所示，該矩形梁的寬度W為0.2 m，高度H為0.5 m，熱傳導(dǎo)系數(shù)k為0.0025，上下邊界溫度為0，熱源項(xiàng)的具體函數(shù)表達(dá)式為：

圖1 矩形柱的幾何模型

HS(y)=2e-y

(21)

相應(yīng)的溫度場(chǎng)的解析解為：

(22)

本算例中將梁的邊界離散為200個(gè)單元，并將利用直線積分法和邊界元法進(jìn)行計(jì)算得到的結(jié)果與解析解進(jìn)行對(duì)比。選取了直線x=0.1 m上的9個(gè)內(nèi)部點(diǎn)作為對(duì)比，具體結(jié)果見表1，結(jié)合圖2，顯而易見，數(shù)值解與解析解高度吻合，進(jìn)而證明了直線積分邊界元法的有效性和精確性。

表1 溫度值對(duì)比

圖2 溫度的解析解與數(shù)值解對(duì)比

3.2 壩體結(jié)構(gòu)的含熱源熱傳導(dǎo)分析

大壩模型具體幾何尺寸如圖3所示，單位為m，熱傳導(dǎo)系數(shù)為1000 W/(m·K)，大壩頂部溫度設(shè)為38.5 ℃，底部溫度為10 ℃，其余邊界均為絕熱狀態(tài)。熱源假設(shè)為常數(shù)0.5 W/m3。

圖3 重力壩幾何模型

為了驗(yàn)證本文方法的正確性，由于缺乏解析解，本文采用有限元法計(jì)算的結(jié)果與本文方法進(jìn)行對(duì)比。如圖4所示，邊界元計(jì)算結(jié)果與有限元計(jì)算結(jié)果十分相近。

圖4 計(jì)算結(jié)果對(duì)比

4 結(jié) 論

本文主要基于邊界元法，進(jìn)行了含熱源的熱傳導(dǎo)研究，由于考慮熱源導(dǎo)致邊界元積分方程中出現(xiàn)了域積分項(xiàng)，為了保證邊界元法的降維優(yōu)點(diǎn)，即只需要對(duì)邊界進(jìn)行離散，本文提出了使用直線積分法進(jìn)行域積分的處理，直線積分法可以基于散度定理，在直角坐標(biāo)系中，將域積分轉(zhuǎn)化為包含直線積分的邊界積分。為了提高直線積分的精度，還可以使用背景網(wǎng)格將直線積分分割為線積分段，之后便可以分別對(duì)直線積分段上的積分點(diǎn)進(jìn)行計(jì)算。結(jié)合數(shù)值驗(yàn)證的第一個(gè)例子可以看出，本文所提出方法的計(jì)算結(jié)果和解析解高度契合，精度很高。此外，還將本文方法運(yùn)用到大壩結(jié)構(gòu)的熱傳導(dǎo)分析當(dāng)中，并與有限元計(jì)算結(jié)果進(jìn)行了對(duì)比，同樣成功地證明了該方法的廣泛適用性。