多網(wǎng)格下解析法計(jì)算石英燈輻射熱流對(duì)比研究

關(guān)鍵詞：多網(wǎng)格； 解析法； 輻射熱流； 有限元法； 石英燈

中圖分類號(hào)：V216.4 文獻(xiàn)標(biāo)識(shí)碼：A DOI：10.19452/j.issn1007-5453.2024.10.008

石英燈加熱裝置是高速飛行器熱結(jié)構(gòu)地面熱強(qiáng)度試驗(yàn)[1]中廣泛采用的加熱裝置。加熱裝置設(shè)計(jì)的合理性直接關(guān)系到結(jié)構(gòu)表面熱載荷施加的準(zhǔn)確性，可以通過(guò)計(jì)算石英燈陣列在結(jié)構(gòu)表面產(chǎn)生的輻射熱流分布情況來(lái)評(píng)估熱載荷的均勻性、峰值水平以及其他分布特征等。

常用的石英燈陣列輻射熱流計(jì)算方法包括解析法、蒙特卡羅法，以及采用基于上述方法的有限元軟件工具等。采用有限元軟件進(jìn)行輻射熱流計(jì)算的缺點(diǎn)主要有：（1）需要對(duì)大量的石英燈元件進(jìn)行實(shí)體建模，需要設(shè)定受熱結(jié)構(gòu)表面與石英燈的復(fù)雜輻射換熱關(guān)系，需要明確受熱結(jié)構(gòu)和石英燈的相關(guān)材料參數(shù)；（2）對(duì)于大規(guī)模石英燈陣列的輻射換熱通常計(jì)算時(shí)間較長(zhǎng)；（3）部分有限元軟件中限制輻射換熱單元總數(shù)，在進(jìn)行大規(guī)模加熱陣列計(jì)算時(shí)，對(duì)石英燈燈絲長(zhǎng)度方向和環(huán)向的網(wǎng)格劃分?jǐn)?shù)量往往只能采取較為粗略的網(wǎng)格劃分，在一定程度上影響計(jì)算結(jié)果的精度?？追步鸬萚2]對(duì)比了區(qū)域法、離散坐標(biāo)法、蒙特卡羅法等幾種經(jīng)典的輻射傳熱計(jì)算方法以及常用的輻射傳熱計(jì)算軟件，討論了幾種輻射熱流場(chǎng)計(jì)算模型的適用性，認(rèn)為可根據(jù)計(jì)算目的和精度要求選擇合適的計(jì)算方法，以實(shí)現(xiàn)對(duì)石英燈陣列熱流場(chǎng)的快速準(zhǔn)確模擬。

蒙特卡羅法通過(guò)光子軌跡追蹤能夠有效模擬光子在石英玻璃管內(nèi)、反射屏作用下的復(fù)雜折、反射過(guò)程，因此計(jì)算精度較高，在輻射熱流計(jì)算中應(yīng)用十分廣泛，國(guó)外很早便建立了采用蒙特卡羅法預(yù)測(cè)石英燈輻射熱流的計(jì)算方法[3-5]，被國(guó)內(nèi)學(xué)者廣泛引用。宋偉浩等[6]利用蒙特卡羅法對(duì)紅外燈陣出射光線進(jìn)行了隨機(jī)抽樣和追跡，建立了受熱平面上光線落點(diǎn)坐標(biāo)數(shù)據(jù)庫(kù)，計(jì)算了受熱平面輻射照度分布情況。夏吝時(shí)等[7]采用蒙特卡羅法對(duì)平板形石英燈加熱器中水冷反光板面積、水冷反光板與燈陣間距離、熱源疏密程度、熱源陣列與材料受熱面間距等因素對(duì)輻射熱場(chǎng)中典型隔熱材料受熱面溫度分布均勻性和熱流密度進(jìn)行了模擬計(jì)算。朱言旦[8-9]基于蒙特卡羅法發(fā)展了石英燈陣輻射熱流分布預(yù)測(cè)方法及計(jì)算程序，以燈陣中各燈功率為優(yōu)化參數(shù)的石英燈陣熱流模擬優(yōu)化設(shè)計(jì)方法，對(duì)某飛行器結(jié)構(gòu)部件迎風(fēng)面氣動(dòng)加熱進(jìn)行了模擬。這些采用蒙特卡羅法的研究工作需要的輸入?yún)?shù)較多、計(jì)算結(jié)果精度直接依賴于樣本數(shù)量大小、計(jì)算時(shí)間相對(duì)較長(zhǎng)等不足。

解析法基于斯忒藩-玻耳茲曼定律、蘭貝特定律等傳熱學(xué)基本定律，求解時(shí)往往需要對(duì)研究對(duì)象進(jìn)行簡(jiǎn)化，使之滿足基本定律的應(yīng)用條件，優(yōu)點(diǎn)是顯式計(jì)算過(guò)程收斂性好，通過(guò)模型簡(jiǎn)化可以實(shí)現(xiàn)較快的計(jì)算速度，同時(shí)可以根據(jù)需求進(jìn)行編程擴(kuò)展，缺點(diǎn)則是不能考慮光線在傳播過(guò)程中復(fù)雜的折、反射現(xiàn)象。宋偉浩等[10]建立了紅外燈陣輻照度分布數(shù)學(xué)模型和受熱平面溫度分布數(shù)學(xué)模型，針對(duì)燈間距離參數(shù)對(duì)受熱平面溫度分布進(jìn)行了仿真與試驗(yàn)。李昕昕等[11]針對(duì)解析法不適用于帶反射屏的石英燈輻射熱流計(jì)算、蒙特卡羅法計(jì)算時(shí)間過(guò)長(zhǎng)的問題，將蒙特卡羅法與解析法相結(jié)合，提出了基于解析法的“等效面光源”求解構(gòu)想，首先采用蒙特卡羅法獲得反射屏作用下的典型熱流場(chǎng)分布結(jié)果，再利用輻射傳熱定律建立關(guān)于面光源各離散點(diǎn)溫度的線性方程組，求解得到等效面光源的溫度分布，之后采用解析法對(duì)等效面光源作用下的熱流場(chǎng)進(jìn)行計(jì)算，實(shí)現(xiàn)了帶反射屏的石英燈加熱裝置的熱流快速計(jì)算。張肖肖等[12]針對(duì)平板形、圓柱形、圓錐形石英燈加熱陣列，提出了基于解析法的熱流快速計(jì)算方法，構(gòu)建了解析傳熱模型，建立了圓柱和圓錐加熱陣列中單根石英燈管照射范圍的判斷依據(jù)，得到了典型燈管排布下熱流分布規(guī)律。這種計(jì)算方法將燈絲半圓柱面作為平面處理，模型得到簡(jiǎn)化，計(jì)算速度顯著提高，但在精度方面有所下降，需要進(jìn)一步探討半圓柱面多網(wǎng)格劃分時(shí)的解析計(jì)算方法。

本文在參考文獻(xiàn)[12]研究工作的基礎(chǔ)上，建立適用于半圓柱面多網(wǎng)格劃分時(shí)的熱流密度解析計(jì)算方法，分析不同網(wǎng)格數(shù)量劃分下熱流分布差異，并與有限元軟件熱流計(jì)算結(jié)果進(jìn)行對(duì)比，以評(píng)估該解析計(jì)算方法的正確性，同時(shí)探討解析法在熱流計(jì)算工程應(yīng)用中的相關(guān)注意事項(xiàng)。

1 多網(wǎng)格時(shí)解析法計(jì)算理論

考慮燈絲面向試驗(yàn)件表面的半圓柱面（見圖1），將燈絲沿長(zhǎng)度方向N等分，對(duì)燈絲半圓柱面按照?qǐng)D2 采用數(shù)量為n的網(wǎng)格進(jìn)行等分，將每一個(gè)弧面視為平面微元處理，把半圓柱面發(fā)出的總熱流平均到每段平面微元上，則每段微元對(duì)外發(fā)出熱流

式中， r 為燈絲半徑；L 為燈絲長(zhǎng)度；N為燈絲沿長(zhǎng)度方向均分?jǐn)?shù)量；σ 為斯忒藩-玻耳茲曼常數(shù)；T為燈絲溫度。

從圖3可以看到，各微元中心位置與接收面上不同位置的接收點(diǎn)之間的相對(duì)關(guān)系呈現(xiàn)出更為復(fù)雜的變化，n=1 時(shí)，可認(rèn)為所有接收點(diǎn)均可受到平面微元的照射，而當(dāng)n>1 時(shí)，對(duì)于非水平位置的平面微元，在接收面上均存在一個(gè)其微元對(duì)應(yīng)半球空間能夠產(chǎn)生照射的極限位置，即微元所在平面與接收面的交點(diǎn)位置（如圖3中的C和C'、D和D'）。

將不同θ位置的各段平面微元對(duì)G點(diǎn)產(chǎn)生的熱流密度累加，便得到了在單根石英燈作用下G點(diǎn)的熱流密度。對(duì)于多根石英燈組成的陣列，可以按照上述過(guò)程逐根計(jì)算每根石英燈對(duì)各接收點(diǎn)產(chǎn)生的熱流密度再累加，便可以得到不同網(wǎng)格數(shù)量劃分下的熱流密度場(chǎng)分布情況。

2 計(jì)算結(jié)果

計(jì)算長(zhǎng)度300mm、直徑2.6mm、距接收面高度50mm的單根石英燈燈絲在300mm×400mm區(qū)域的熱流分布情況，接收面長(zhǎng)度和寬度方向每10mm選取一個(gè)節(jié)點(diǎn)，燈絲沿長(zhǎng)度方向等分為100 段，燈絲溫度取2000K，計(jì)算得到不同網(wǎng)格數(shù)量下的熱流分布云圖如圖5 所示，圖中白色虛線為燈絲在接收面上的投影。可以看到當(dāng)n 從1 增加到4 時(shí)，熱流分布以及熱流峰值水平呈現(xiàn)出相對(duì)顯著的變化，而當(dāng)n 從4繼續(xù)增加后，熱流分布以及熱流峰值水平相對(duì)變化不明顯。

提取上述云圖中y=200 以及x=150 兩條直線上各節(jié)點(diǎn)以及坐標(biāo)為（150，200）的中心點(diǎn)的熱流密度計(jì)算結(jié)果，如圖6 所示?？梢钥吹剑S著n 的增大，熱流密度峰值逐漸減小，當(dāng)n 從1 增加到2 時(shí)，熱流密度峰值變化尤為明顯，當(dāng)n 從4繼續(xù)增加后，峰值水平逐漸趨于穩(wěn)定；與中心點(diǎn)處相反，坐標(biāo)為（0，200）和（300，200）的左右邊緣處熱流密度，當(dāng)n=1 時(shí)取得最小值；而坐標(biāo)為（150，0）和（150，400）的上下邊緣處熱流密度隨著n 的增大，呈現(xiàn)出與中心點(diǎn)處相似的變化趨勢(shì)。這是由于當(dāng)n=1 時(shí)燈絲所有輻射能量集中于一個(gè)矩形窄帶內(nèi)，根據(jù)蘭貝特（Lambert）定律，此時(shí)燈絲對(duì)正下方區(qū)域的定向輻射強(qiáng)度達(dá)到最高值，而對(duì)左右邊緣處的輻射作用受較小的可見面積影響相對(duì)最弱；而隨著n 的增加，燈絲的輻射能量均勻分散于沿環(huán)向的各矩形窄帶內(nèi)，燈絲對(duì)正下方區(qū)域的輻射強(qiáng)度逐漸降低，但由于兩者間距離最短，所以正下方區(qū)域的熱流密度仍然達(dá)到當(dāng)前分布的最高值，而對(duì)左右邊緣處的輻射作用受到環(huán)向窄帶對(duì)其可見面積增大的影響相對(duì)n=1 有了一定增強(qiáng)，但受到可對(duì)其照射的窄帶數(shù)量降低的影響，增強(qiáng)幅度較??；隨著n 的增加，環(huán)向窄帶對(duì)上下邊緣處的可見面積相對(duì)n=1 進(jìn)一步降低，因此上下邊緣處受到的輻射作用當(dāng)n=1時(shí)達(dá)到最強(qiáng)。

3 有限元計(jì)算結(jié)果對(duì)比

3.1 熱流計(jì)算結(jié)果對(duì)比

目前，商用有限元軟件在石英燈輻射熱流計(jì)算中被廣泛采用，將解析法計(jì)算結(jié)果與相同參數(shù)條件下的有限元計(jì)算結(jié)果進(jìn)行對(duì)比，通過(guò)結(jié)果的一致性可以驗(yàn)證解析法的正確性。采用與第2 節(jié)中相同尺寸、高度和溫度的燈絲以及相同大小的接收面，接收面網(wǎng)格尺寸為10mm×10mm，燈絲沿長(zhǎng)度方向等分為100個(gè)網(wǎng)格，按照面面輻射法在有限元軟件計(jì)算不同環(huán)向網(wǎng)格數(shù)量下的熱流密度分布，計(jì)算模型如圖7 所示。同樣提取橫豎中軸線上的熱流分布結(jié)果，與解析法計(jì)算結(jié)果進(jìn)行對(duì)比，如圖8 所示。從圖8 中可以看到，當(dāng)n=1 時(shí)（此時(shí)燈絲有限元模型與解析法模型均為矩形窄帶），解析法與有限元軟件計(jì)算結(jié)果呈現(xiàn)出顯著差異，峰值熱流水平相對(duì)誤差達(dá)到61.7%；隨著n 的增加（此時(shí)燈絲有限元模型與解析法模型均為半圓柱面，并且環(huán)向網(wǎng)格數(shù)相同），兩種方法的熱流計(jì)算結(jié)果的差異逐漸縮小，當(dāng)n>4時(shí)，計(jì)算結(jié)果一致性相對(duì)保持穩(wěn)定；與解析法計(jì)算結(jié)果不同，有限元軟件計(jì)算得到的中心點(diǎn)處熱流密度隨著n 的增大呈現(xiàn)出小幅的增大。

因此，燈絲對(duì)外輻射總能量差異應(yīng)是解析法與有限元熱流計(jì)算結(jié)果差異的最主要原因。α 隨n 的變化如圖9 所示，從圖9 可以看到，隨著n 的增大，輻射總能量之比逐漸趨近于1；當(dāng)n=1 時(shí)，解析法中輻射總能量是有限元軟件中的1.57 倍，在相同燈絲溫度情況下，有限元軟件中采用燈絲直徑寬度的矩形窄帶來(lái)計(jì)算石英燈輻射傳熱會(huì)顯著降低結(jié)構(gòu)溫度響應(yīng)的精度；當(dāng)n=3 時(shí)，輻射總能量之比達(dá)到1.047，n=4 時(shí)，輻射總能量之比達(dá)到1.026，因此建議在有限元建模時(shí)對(duì)燈絲半圓柱至少采取3～4個(gè)的網(wǎng)格進(jìn)行劃分。

為驗(yàn)證解析法的有效性，將有限元計(jì)算中的燈絲溫度結(jié)合式（9）進(jìn)行調(diào)整，以保證二者在計(jì)算過(guò)程中燈絲對(duì)外輻射總能量一致，進(jìn)而通過(guò)對(duì)比驗(yàn)證解析法計(jì)算結(jié)果的正確性。將有限元計(jì)算中的燈絲溫度T調(diào)整為式（10）中的T

式中，T為調(diào)整后的有限元計(jì)算中燈絲溫度；T 為給定的石英燈燈絲溫度。

3.3 調(diào)整燈絲溫度后的熱流計(jì)算結(jié)果對(duì)比

對(duì)有限元模型中的燈絲溫度按照式（10）調(diào)整后，重新計(jì)算并對(duì)比解析法與有限元計(jì)算熱流密度結(jié)果，如圖10所示。從圖10可以看到，經(jīng)過(guò)燈絲溫度調(diào)整、保證輻射總能量一致后，有限元計(jì)算結(jié)果與解析法計(jì)算結(jié)果的一致性相對(duì)調(diào)整前得到顯著改善，n=1 時(shí)中心點(diǎn)熱流密度相對(duì)誤差僅有3.23%。其他產(chǎn)生誤差的原因可能是有限元軟件中熱流密度為單元計(jì)算結(jié)果，此處的節(jié)點(diǎn)熱流密度結(jié)果是將其相鄰單元的結(jié)果進(jìn)行平均得到，而解析法計(jì)算結(jié)果則是該節(jié)點(diǎn)位置的理論計(jì)算結(jié)果，兩者存在一定偏差。

4 結(jié)論

針對(duì)石英燈輻射熱流計(jì)算，本文提出了適用于半圓柱面多網(wǎng)格劃分時(shí)的解析計(jì)算方法，對(duì)比了與有限元軟件熱流計(jì)算結(jié)果的差異，確定了差異產(chǎn)生原因并建立了調(diào)整方法，驗(yàn)證了解析計(jì)算方法的正確性，得到了以下結(jié)論：

（1）采用解析法計(jì)算石英燈熱流密度分布時(shí)，采用不同的網(wǎng)格數(shù)量劃分燈絲半圓柱面會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生影響，當(dāng)n≥4后，熱流密度峰值水平逐漸趨于穩(wěn)定。

（2）n=1 時(shí)采用解析法計(jì)算石英燈熱流密度雖然計(jì)算精度相對(duì)較差，但由于燈絲網(wǎng)格數(shù)量較少且不存在照射范圍判斷問題，在同等計(jì)算條件下能夠獲得最快的計(jì)算速度，因此在大規(guī)模石英燈陣列的輻射熱流預(yù)估中仍然有其工程應(yīng)用價(jià)值；當(dāng)n 達(dá)到3 以上時(shí)，解析法計(jì)算結(jié)果可以滿足工程中對(duì)石英燈熱流仿真計(jì)算精度的要求。

（3）在相同燈絲溫度情況下，有限元軟件中采用燈絲直徑寬度的矩形窄帶來(lái)計(jì)算石英燈輻射傳熱會(huì)顯著降低結(jié)構(gòu)溫度響應(yīng)的精度，因此在有限元建模時(shí)對(duì)燈絲半圓柱面至少采取3～4個(gè)的網(wǎng)格進(jìn)行劃分。

（4）燈絲對(duì)外輻射總能量差異是解析法與有限元熱流計(jì)算結(jié)果差異的最主要原因，將燈絲溫度按照總能量比進(jìn)行調(diào)整后，解析法與有限元計(jì)算結(jié)果的一致性得到大幅改善，證明了本文提出的多網(wǎng)格時(shí)熱流密度解析計(jì)算方法的正確性。

（5）采用本文提出的解析法計(jì)算的熱流密度值可作為熱流載荷邊界施加于有限元分析模型中，從而省去在有限元軟件中對(duì)石英燈建模的過(guò)程，同時(shí)可以不受部分有限元軟件對(duì)輻射換熱單元總數(shù)的限制，能夠滿足大尺度結(jié)構(gòu)與大規(guī)模石英燈陣列輻射傳熱下結(jié)構(gòu)溫度響應(yīng)快速計(jì)算的需求。