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

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

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

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

而當面臨考慮熱源的熱傳導問題時，在邊界積分方程中，會出現域積分，如果不加以處理，則需要進行單元離散，使得邊界元喪失了降維的優(yōu)勢。因此，尋找一種能夠將域積分轉化為邊界積分的方法是十分必要的。

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

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

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

結合前人研究，可知控制方程為：

(1)

式中：k為熱傳導系數；B(x)為熱源項；θ為溫度；x為研究域Ω內的點，其組成元素為xi。

邊界情況可設為：

(2)

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

(3)

(4)

(5)

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

r=‖x-y‖

(6)

可以看到，式中含有一項域積分：

(7)

需要利用直線積分法將域積分轉化為邊界積分。

2 直線積分法

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

(8)

假設函數f(x,y)可以從矢量場的散度得到，具體表達為：

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

(9)

式中：div為散度運算符；F為向量場，可表達為：

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

(10)

式中：ei為直角坐標系中的一組基向量。

通過運用散度定理，可以將域積分轉化為以下形式：

(11)

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

(12)

或

(13)

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

(14)

(15)

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

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

(16)

(17)

式(17)的計算一般都可以采用數值計算的方式，現在可以將計算域用背景網格覆蓋并逐級等分為小網格，進而可以將積分線截斷成為一段段的積分線，則式(17)可以改寫為：

(18)

式中：Lk為背景網格劃分后的直線積分線段。

引入背景網格之后，利用每個網格內的積分線段上的高斯積分點就可以直接進行計算并且可以獲得很高的精度。進一步地，域積分可轉化為：

(19)

最后，將直線積分法應用于邊界積分方程中的域積分的計算，具體轉換后的表達式為：

(20)

3 數值驗證

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

3.1 矩形柱的含熱源熱傳導分析

如圖1所示，該矩形梁的寬度W為0.2 m，高度H為0.5 m，熱傳導系數k為0.0025，上下邊界溫度為0，熱源項的具體函數表達式為：

圖1 矩形柱的幾何模型

HS(y)=2e-y

(21)

相應的溫度場的解析解為：

(22)

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

表1 溫度值對比

圖2 溫度的解析解與數值解對比

3.2 壩體結構的含熱源熱傳導分析

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

圖3 重力壩幾何模型

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

圖4 計算結果對比

4 結 論

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