關鍵詞：瞬態，熱傳導，有限元求解器，三角形單元

熱傳遞有三種方式：熱傳導、熱對流、熱輻射。就熱傳導問題而言，無論是結構力學還是流體力學都會涉及，兩邊都沒拿它當外人。

前面的文章提到過，結構力學的有限元發展地非常成熟，大部分的剛度矩陣在文獻里面都推導好了。而流體力學的很多單元類型的有限元方程，可能需要自行推導完成。在熱傳導問題中，我采用加權余量法進行處理，推導出了符合結構力學有限元文獻中給出的剛度矩陣，殊途同歸。

實際上，傳統的結構力學有限元三大控制方程：幾何方程、物理方程、平衡方程。幾何方程描述位移-應變關系，物理方程描述應力-應變關系，平衡方程描述內應力-外載荷關系。傳熱問題從控制方程角度，更偏向流體力學（能量方程）。但是熱對于結構變形太重要了，因此結構有限元必須要把傳熱問題解決掉。

從結構力學跨到流體力學，在有限元方法中，流體力學控制方程左邊的矩陣都可以用剛度矩陣去看待它。控制方程的右邊的列陣，都可以用載荷的角度去看待，對于第二類邊界條件，則可以分成左側矩陣的修正+右側列陣的載荷組合。有些文獻上，用所謂的“內部單元方程”、“邊界單元方程”的描述，會增加我們的困惑，可以不必糾結在此。

控制方程

二維瞬態熱傳導控制方程如下：

這個方程里面的常數有密度、比熱容、導熱系數。

三種邊界條件：

（1） 已知邊界溫度值，屬于第一類邊界條件，它的處理就和結構有限元里面的位移以一樣，可以用置大數法對方程左邊的矩陣進行約束處理。

（2） 已知邊界熱流密度，屬于第二類邊界條件，作為熱源。可以類比到結構有限元里面的均布載荷。

（2） 已知邊界對流換熱系數和接觸環境溫度，也屬于第二類邊界條件。這個邊界條件在處理的時候，需要進行拆分，一部分放到左側單元矩陣，一部分作為右側的載荷。

有限元思路

這部分在結構有限元教材中介紹的比較多，流程：

（1） 根據單元類型，確定插值函數。此時單元溫度用權函數表達。

（2） 采用伽遼金方法，權函數=插值函數，控制方程與權函數相乘，積分取0。

（3） 在每個單元域內，方程轉換為權函數的積分形式，最終形成單元矩陣。

單元方程

使用三角形線性單元對應的插值函數：

有些教材中，會把面積項提取出來，寫成以下這種形式，所以有的教材上剛度矩陣結果用a、b、c表達的時候，會存在差異，但是本質都是一樣的。

最終單元方程如下，其中M是熱容矩陣，K是傳導矩陣，F是熱載荷。

熱容矩陣乘的是溫度的導數。在瞬態問題的求解中，導數項可以寫成前后時間變量差值與時間間隔的比值：

代入后得到如下形式：

求解思路

在求解過程中，把Tn+1當作未知量，Tn作為已知量。這樣在每個時間點，求解方法和結構有限元方法一致。

初始時候，可以指定一個溫度作為全域已知初始溫度，然后在迭代過程中，Tn和Tn+1會逐漸接近，達到收斂狀態。

案例效果

設計案例如下，同時包含對流換熱邊界條件和熱流，時間總長10000s，每步時間間隔50s。

自研求解器和商用軟件結果對比如下，從結果可以看出，自研求解器結果與商用軟件結果一致。

自研求解器結果：最終溫度分布

商用軟件結果：最終溫度分布

自研求解器結果：平均溫度時間曲線

商用軟件結果：平均溫度時間曲線