高效偉 丁金興 劉華雩

(大連理工大學(xué)航空航天學(xué)院,工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,大連 116024)

本文創(chuàng)建一種全新的數(shù)值計(jì)算方法—有限線法,并將其用于求解流體-固體一體化耦合傳熱分析.常用的有限元法是基于體單元的離散方法,有限容積法是在控制體面上運(yùn)算的方法,邊界元法是在邊界面上離散的方法,無(wú)網(wǎng)格法等是由計(jì)算點(diǎn)周圍的散點(diǎn)構(gòu)建計(jì)算格式的方法.本文提出的有限線法是一種基于有限條線段離散的方法,在每個(gè)點(diǎn)只需要兩條或三條直線或曲線構(gòu)成的線系,則可建立任意高階的算法格式.創(chuàng)新性思想是: 通過(guò)采用沿線段求方向?qū)?shù)的技術(shù),由一維拉格朗日插值公式,建立二維和三維問(wèn)題的任意高階線系的空間導(dǎo)數(shù),并以此直接由問(wèn)題的控制偏微分方程與邊界條件建立離散的系統(tǒng)方程組.有限線法理論簡(jiǎn)單、通用性強(qiáng),能以統(tǒng)一的格式求解固體與流體力學(xué)等偏微分方程邊值問(wèn)題.在流體方程中,擴(kuò)散項(xiàng)采用中心線系保證高精度計(jì)算,而對(duì)流項(xiàng)則采用迎風(fēng)線系體現(xiàn)其方向性特征.本文將給出有限線法求解流固耦合傳熱問(wèn)題的幾個(gè)算例分析,驗(yàn)證其正確性與有效性.

1 引言

根據(jù)算法離散與操作維度的不同,常用的數(shù)值計(jì)算方法可歸納為4 類: 基于體單元的計(jì)算方法(如有限元法[1,2]、有限塊法[3,4]、單元微分法[5,6]等),基于面單元的方法(如有限容積法[7,8]、邊界元法[9,10]、邊界面法[11]等),基于線單元的方法(如有限元線法[12]、交叉線法[13]等),以及基于散點(diǎn)計(jì)算的方法(如無(wú)網(wǎng)格法[14-16]、基本解法[17]、奇異邊界點(diǎn)法[18]等).

在基于體單元的計(jì)算方法中最具代表性的方法是有限元法(FEM).用其分析問(wèn)題時(shí)將計(jì)算域劃分成有限個(gè)體單元,物理量的表征、微分是通過(guò)與計(jì)算域維度一樣的體單元實(shí)現(xiàn)的.鑒于其強(qiáng)大的解題能力,FEM 已被用于解決各種工程問(wèn)題,特別是固體力學(xué)問(wèn)題,這主要?dú)w因于其使用了具有良好性能的等參單元進(jìn)行幾何離散和物理量插值[19],使得計(jì)算結(jié)果非常穩(wěn)定.此外,FEM 所需要的單元和節(jié)點(diǎn)數(shù)遠(yuǎn)少于有限容積法(FVM)和有限差分法(FDM).FEM的弱點(diǎn)主要體現(xiàn)在: 1)需使用變分原理或虛功原理來(lái)建立計(jì)算格式,對(duì)于不同的問(wèn)題所建立的表達(dá)式不同,從而不便于建立多場(chǎng)耦合問(wèn)題的統(tǒng)一算法;2)用傳統(tǒng)的FEM 求解含有對(duì)流項(xiàng)的問(wèn)題時(shí),不易實(shí)施迎風(fēng)格式與激波捕捉計(jì)算.為此,近年來(lái)間斷伽遼金有限元法(DGFEM)在流體力學(xué)中得到了快速發(fā)展[20],其在每個(gè)單元內(nèi)用獨(dú)立定義的節(jié)點(diǎn)進(jìn)行高階函數(shù)逼近,因而可以較容易實(shí)現(xiàn)迎風(fēng)計(jì)算,然而DGFEM 計(jì)算效率很低,此外如果單元內(nèi)有激波的存在,內(nèi)部間斷使得收斂性和魯棒性大幅下降.

面單元計(jì)算方法是指主要在面上進(jìn)行運(yùn)算的方法.將FVM 歸入面單元計(jì)算方法在維度上似乎不太協(xié)調(diào),但是科學(xué)的.這是因?yàn)镕VM 是一種將問(wèn)題的控制微分方程的守恒形式(全微分形式)使用高斯散度公式將控制體積分轉(zhuǎn)換成在控制體面上進(jìn)行積分運(yùn)算的算法[7].FVM 是一種通量守恒性最好、激波間斷面最容易捕捉的數(shù)值方法,也由此成了目前計(jì)算流體力學(xué)中的主流方法.BEM 是一種典型的面單元計(jì)算法,是通過(guò)兩次分部積分和使用高斯散度公式,并取權(quán)函數(shù)為問(wèn)題基本解而形成邊界積分方程的算法.BEM 幾何離散簡(jiǎn)單,能很好地模擬斷裂力學(xué)中的應(yīng)力集中現(xiàn)象,并能有效求解無(wú)限域與半無(wú)限域問(wèn)題[21].BEM的弱點(diǎn)是在求解有源(熱源、體積力、化學(xué)反應(yīng)生成項(xiàng)等)、非線性、非均質(zhì)問(wèn)題時(shí),會(huì)有域積分出現(xiàn)在積分方程中,削弱了其只在邊界上離散的優(yōu)點(diǎn).

基于線單元的計(jì)算方法籠統(tǒng)地講可以包括有限元線法[12]以及梁?jiǎn)卧?、桿單元[22]等的算法.前者是借助于能量泛函的變分,將控制偏微分方程半離散化為用結(jié)線表示的常微分方程組求解法,而后者則是通過(guò)一維離散模擬具有一維維度的構(gòu)件的方法.真正通過(guò)線單元的離散實(shí)現(xiàn)高維(二維和三維)問(wèn)題的微分方程的求解方法是作者最近提出的交叉線法(CLM)[13].該方法通過(guò)在每個(gè)點(diǎn)建立由兩條(二維問(wèn)題)或三條(三維問(wèn)題)相交的線段,并通過(guò)形函數(shù)構(gòu)造法建立多維問(wèn)題的形函數(shù),然后采用微分法則,導(dǎo)出各階偏導(dǎo)數(shù)計(jì)算式,最后將其代入控制微分方程與邊界條件,建立問(wèn)題的離散系統(tǒng)方程組.CLM 屬于配點(diǎn)型的強(qiáng)形式計(jì)算方法,目前是通過(guò)單元的形式在自由單元法中使用[13],用到的線系可以是直線,也可以是曲線.因而,該方法具有幾何離散簡(jiǎn)單、適應(yīng)性強(qiáng)、編程方便、系數(shù)矩陣稀疏等特點(diǎn),是一種具有開(kāi)發(fā)潛力的計(jì)算方法.然而,根據(jù)單元形函數(shù)構(gòu)造法(即滿足Delta符號(hào)性質(zhì)的解[13])很難構(gòu)造三階及其以上的形函數(shù),這給流體力學(xué)中經(jīng)常用到高階格式的算法帶來(lái)了限制.

基于散點(diǎn)的計(jì)算方法有無(wú)網(wǎng)格法(MLM)等[14-18],該類方法不需要有規(guī)則節(jié)點(diǎn)分布的單元,只需一些任意分布的散點(diǎn)則可建立解題算法.MLM的形函數(shù)及其偏導(dǎo)數(shù)是通過(guò)最小二乘法等技術(shù)在一定范圍(緊支域)內(nèi)的點(diǎn)集構(gòu)建的.基本解法與奇異邊界點(diǎn)法等則是通過(guò)問(wèn)題微分算子的基本解在虛擬邊界或真實(shí)邊界節(jié)點(diǎn)上配置方程進(jìn)行求解的.MLM在幾何適應(yīng)方面具有獨(dú)特的優(yōu)勢(shì),因而在復(fù)雜幾何模擬與動(dòng)邊界問(wèn)題方面得到了廣泛的應(yīng)用.MLM雖然有使用靈活、幾何適應(yīng)性強(qiáng)的優(yōu)點(diǎn),但也有計(jì)算效率低、結(jié)果穩(wěn)定性差以及不易施加邊界條件與各向異性物性的弱點(diǎn).

值得介紹的一種新的配點(diǎn)型方法是作者近年來(lái)提出的自由單元法(FrEM)[23,24].該方法吸收了FEM 與配點(diǎn)型MLM的優(yōu)點(diǎn),即: 在函數(shù)表征方面,采用FEM 中有規(guī)則節(jié)點(diǎn)分布的等參單元,因而具有幾何適應(yīng)性強(qiáng)、計(jì)算效率高、結(jié)果穩(wěn)定性好、容易施加各向異性物性以及邊界條件的優(yōu)點(diǎn);在算法方面,采用配點(diǎn)法技術(shù),直接由控制方程與邊界條件逐點(diǎn)產(chǎn)生系統(tǒng)方程.FrEM 在每個(gè)點(diǎn)只需要一個(gè)和周圍自由選擇的節(jié)點(diǎn)形成的獨(dú)立的等參單元,因而不需要考慮物理量在單元之間的相互連接關(guān)系以及導(dǎo)數(shù)連續(xù)性問(wèn)題.此外,單元的自由性也使其更容易構(gòu)造高階迎風(fēng)格式計(jì)算流體力學(xué)方程中的對(duì)流項(xiàng).FrEM 已在固體力學(xué)[25]、流體力學(xué)[24]與傳熱學(xué)[26]中得到了成功的應(yīng)用.然而,由于高階自由單元是由結(jié)構(gòu)化節(jié)點(diǎn)形成的[23],因此三維問(wèn)題的高階單元需要較多的節(jié)點(diǎn),這給幾何建模與數(shù)據(jù)存儲(chǔ)帶來(lái)了很大的不便.

本文在作者多年來(lái)研究BEM、FEM、配點(diǎn)型MLM 等方法的基礎(chǔ)上,提出了一種全新的數(shù)值方法—有限線法(FLM)[27].該方法在每個(gè)點(diǎn)只需要1—3 條交叉線段(稱為線系)則可建立一維到三維問(wèn)題的算法格式.由于每條線段用拉格朗日多項(xiàng)式插值公式構(gòu)建函數(shù)表征關(guān)系,因而可以很方便地建立任意高階的算法.FLM 具有CLM 和FrEM的全部?jī)?yōu)點(diǎn),還有類似于FDM 容易使用的性能,并有高階線系所用節(jié)點(diǎn)少、容易施加不同物性與邊界條件、以及容易構(gòu)建迎風(fēng)格式的特點(diǎn),因此既適合于求解流體力學(xué)問(wèn)題、也適合求解固體力學(xué)問(wèn)題.

一些現(xiàn)代裝備結(jié)構(gòu),如發(fā)動(dòng)機(jī)渦輪等,一直處于高溫服役環(huán)境,因而使用冷卻技術(shù)提高其耐高溫能力是經(jīng)常使用的重要手段,結(jié)構(gòu)內(nèi)部設(shè)置冷卻通道是最常用的熱防護(hù)形式,對(duì)其進(jìn)行流體域-固體域間耦合傳熱分析是結(jié)構(gòu)設(shè)計(jì)的重要依據(jù).本文將系統(tǒng)介紹FLM 在求解冷卻通道結(jié)構(gòu)流固域間耦合傳熱方面的應(yīng)用,為現(xiàn)代裝備結(jié)構(gòu)設(shè)計(jì)提供有力的分析手段.

2 流體與固體中的傳熱偏微分方程邊值問(wèn)題

流體與固體中的熱傳導(dǎo)微分方程可統(tǒng)一表示為[5,24]

式中,ρ為流體密度,cp為定壓比熱,vi為流速,T為溫度,λij為導(dǎo)熱率,s為熱源;式中的重復(fù)指標(biāo)表示求和,Ω為計(jì)算域,對(duì)于固體vi為0.(1) 式中的左端為對(duì)流項(xiàng),右端第一項(xiàng)為擴(kuò)散項(xiàng)[7].

相關(guān)的邊界條件可表示為

其中,?Ω為Ω的邊界,ni為邊界外法線方向矢量,h為換熱系數(shù),T∞為環(huán)境溫度.

從 (1)—(4) 式可以看出,控制微分方程與邊界條件涉及到的關(guān)鍵項(xiàng)是溫度對(duì)空間總體坐標(biāo)的一階與二階偏導(dǎo)數(shù),本文將用有限線法建立各階偏導(dǎo)數(shù)計(jì)算式.

3 計(jì)算空間偏導(dǎo)數(shù)的有限線法

作者在文獻(xiàn)[5]中導(dǎo)出了一般拉格朗日等參單元形函數(shù)對(duì)總體坐標(biāo)的偏導(dǎo)數(shù)解析表達(dá)式,但可以證明其一階偏導(dǎo)數(shù)只與通過(guò)計(jì)算點(diǎn)的交叉線上的節(jié)點(diǎn)有關(guān)[13],因此本文采用交叉線系建立各階空間偏導(dǎo)數(shù)計(jì)算式.

3.1 交叉線系的建立

本文提出的FLM 是一種配點(diǎn)型的計(jì)算方法,計(jì)算時(shí)與FDM 和MLM 類似,首先通過(guò)網(wǎng)格生成技術(shù)[28],將計(jì)算域劃分成一系列配置點(diǎn)(CP),圖1給出了一個(gè)二維問(wèn)題的節(jié)點(diǎn)劃分示例.

圖1 將計(jì)算域劃分成一系列配置點(diǎn)Fig.1.Discretizing computational domain into a series of collocation points.

然后,對(duì)每個(gè)CP,建立由兩條(二維問(wèn)題)或三條(三維問(wèn)題)線段組成的交叉線系.圖2 和圖3分別展示了一個(gè)內(nèi)部點(diǎn)的二維(d=2)和三維(d=3)的交叉線系.原理上,線系中的線段相互間越垂直計(jì)算精度越高,但只要不共線(二維線系)和不共面(三維線系)則可獲得滿意的結(jié)果.

對(duì)于位于邊界和角點(diǎn)上的CP,則需要將圖2和圖3 中的交叉點(diǎn)移到邊界和角點(diǎn)上,具體做法可參見(jiàn)文獻(xiàn)[13].

圖2 兩條線組成的二維(d=2)線系Fig.2.Line set of 2D (d=2) consisting of two lines.

圖3 三條線組成的三維(d=3)線系Fig.3.Line set of 3D (d=3) consisting of three lines.

3.2 基于線系的空間導(dǎo)數(shù)計(jì)算公式

在線系的每一條線上,定義m個(gè)節(jié)點(diǎn)(如圖2和圖3 所示),并由下列拉格朗日多項(xiàng)式插值公式表達(dá)坐標(biāo)x與物理量u的函數(shù)變化關(guān)系[5,23]:

式中,l為從節(jié)點(diǎn)1 開(kāi)始計(jì)算的線段長(zhǎng)度,uα為第α個(gè)節(jié)點(diǎn)的u值.Lα為拉格朗日多項(xiàng)式插值函數(shù):

其中,lα為第α個(gè)節(jié)點(diǎn)的長(zhǎng)度.以二維線系為例,標(biāo)記函數(shù)u對(duì)總體坐標(biāo)xi的一階偏導(dǎo)數(shù)分量為?u/?x1和?u/?x2,則沿線段l的方向?qū)?shù)?u/?l可表示為

對(duì)于每一條線可建立一個(gè)上述方程,當(dāng)d條線考慮后可得到如下的矩陣方程:

對(duì)(9)式求逆,并考慮到(6)式,則可得如下的分量表達(dá)式:

為了方便使用,可將(11)式簡(jiǎn)寫(xiě)為

式中,重復(fù)指標(biāo)α表示對(duì)M1個(gè)節(jié)點(diǎn)求和,xc表示線系中各線交叉處(即配點(diǎn))的坐標(biāo)稱為一階導(dǎo)數(shù)算子.需要說(shuō)明的是,矩陣J是一個(gè)2×2 或3×3的矩陣,可以很容易求得其逆矩陣的解析表達(dá)式.

對(duì)于二階導(dǎo)數(shù),可以在一階導(dǎo)數(shù)公式上再次作用導(dǎo)數(shù)算子,得到:

需要說(shuō)明的是,(13)式與(15)式中重復(fù)指標(biāo)α求和的節(jié)點(diǎn)數(shù)量是不一樣的.圖4 展示了一個(gè)配置點(diǎn)(節(jié)點(diǎn)3)采用5×3 線系離散中,一階導(dǎo)數(shù)相關(guān)的7 個(gè)節(jié)點(diǎn)(由實(shí)心點(diǎn)標(biāo)記的節(jié)點(diǎn))與21 個(gè)二階導(dǎo)數(shù)相關(guān)的節(jié)點(diǎn)(由實(shí)心點(diǎn)+空心圓點(diǎn)標(biāo)記的節(jié)點(diǎn))的示意圖.對(duì)于更高階的導(dǎo)數(shù)算子可由類似的遞歸方法形成[27].

圖4 一階與二階導(dǎo)數(shù)相關(guān)的節(jié)點(diǎn)Fig.4.Related nodes of the 1st and 2nd order partial derivatives.

4 流固耦合傳熱問(wèn)題中的有限線法公式

將空間導(dǎo)數(shù)計(jì)算公式(13)和(15)直接代入流固耦合傳熱方程(1)及其邊界條件(2)—(4),則可建立耦合問(wèn)題的代數(shù)離散方程.不過(guò),流體傳熱方程中存在對(duì)流項(xiàng),需要用迎風(fēng)格式才能精確求解,而擴(kuò)散項(xiàng)用中心格式才能得到穩(wěn)定的計(jì)算結(jié)果.本研究通過(guò)采用兩套線系的技術(shù)解決此問(wèn)題,即使用迎風(fēng)線系體現(xiàn)計(jì)算點(diǎn)上游對(duì)其的影響,用中心線系保證擴(kuò)散項(xiàng)的穩(wěn)定性與精度.圖5 展示了采用5×4 節(jié)點(diǎn)線段形成的兩套線系的節(jié)點(diǎn)分布,其中紅色的節(jié)點(diǎn)屬于中心線系(central line set),黑色的屬于迎風(fēng)線系(upwind line set).

圖5 同一配置點(diǎn)的迎風(fēng)與中心線系Fig.5.Upwind and central line sets at a same collocation point.

采用 (13) 式和(15) 式可將方程(1)—(4)離散成下列形式:

式中分別為基于迎風(fēng)線系和中心線系計(jì)算的導(dǎo)數(shù)算子,其中的q為

耦合計(jì)算時(shí),將流體域與固體域進(jìn)行一體化離散.對(duì)于固體域,(17)式的第1 個(gè)方程中的第1 項(xiàng)(對(duì)流項(xiàng))為0.

將上述方程中的配置點(diǎn)xc作用于流體域與固體域的所有配置點(diǎn),并在流-固界面上使用下列溫度協(xié)調(diào)條件與通量平衡方程(其中上標(biāo)f 與s 分別代表流體與固體):

則可形成如下的一體化分析矩陣方程:

求解上述方程則可獲得所有節(jié)點(diǎn)的溫度值.需要說(shuō)明的是,由于每個(gè)配置點(diǎn)的線系所涉及到的節(jié)點(diǎn)數(shù)比傳統(tǒng)的數(shù)值方法FEM 和FVM 少得多,因此 (20) 式中的系數(shù)矩陣A是一個(gè)非常稀疏的非對(duì)稱矩陣.

5 算例分析

根據(jù)本文所推導(dǎo)的公式,采用FORTRAN 語(yǔ)言編制了有限線法程序,程序中采用了分區(qū)技術(shù)[28],可對(duì)復(fù)雜問(wèn)題進(jìn)行計(jì)算分析.下面對(duì)現(xiàn)代高溫裝備中經(jīng)常使用的冷卻管道結(jié)構(gòu)中的一些流固域間耦合傳熱問(wèn)題進(jìn)行分析,用以驗(yàn)證本文所述方法的正確性與有效性.計(jì)算中,每個(gè)配置點(diǎn)的每條迎風(fēng)線用二階精度,中心線用四階精度.

5.1 單根管道結(jié)構(gòu)流動(dòng)傳熱分析

本文的第一個(gè)算例是用于驗(yàn)證本文提出的FLM的正確性與有效性.算例的幾何尺寸與邊界條件如圖6 所示,管道的上下外表面處于環(huán)境溫度為900 K的對(duì)流換熱環(huán)境中,流固界面的換熱系數(shù)為1100 W·(m2·K)—1,流體域的物性為:λ=0.6 W·(m·K)—1,ρ=1000 kg/m3,cp=4200 J·(kg·K)—1;固體域?yàn)棣?200 W·(m·K)—1.

圖6 單根管道結(jié)構(gòu)尺寸與邊界條件Fig.6.Dimensions and B.C.of a single channel structure.

為了提供參考結(jié)果進(jìn)行比較,本文用FLUENT軟件采用較為稠密的網(wǎng)格對(duì)該問(wèn)題進(jìn)行計(jì)算,其中在x方向(水平方向)和固體域,網(wǎng)格為等間距分布,流體域靠近壁面網(wǎng)格逐漸加密.在x方向劃分了151 個(gè)點(diǎn),在垂直方向: 固體域20 個(gè)點(diǎn),流體域40 個(gè)點(diǎn).

5.1.1 FLM 使用均勻網(wǎng)格的網(wǎng)格收斂性分析

首先基于均勻網(wǎng)格對(duì)該問(wèn)題進(jìn)行FLM 計(jì)算分析,考察其網(wǎng)格收斂性.計(jì)算發(fā)現(xiàn),結(jié)果只對(duì)流體域中的垂直方向的網(wǎng)格敏感.因此,在FLM 計(jì)算中,在x方向流固域都使用51 個(gè)點(diǎn);在垂直方向,固體域使用7 個(gè)點(diǎn),流體域分別使用21、41 和71 個(gè)點(diǎn)進(jìn)行計(jì)算,這3 種情況分別標(biāo)記為FLMfine1,FLM-fine2 和 FLM-fine3.圖7 和圖8 分別給出了流速為v=1 m/s 時(shí)3 種網(wǎng)格的計(jì)算結(jié)果.可以看出,網(wǎng)格FLM-fine3 收斂到了FLUENT的計(jì)算結(jié)果.此外,變化趨勢(shì)顯示,FLM 具有很好的網(wǎng)格收斂性.

圖8 不同網(wǎng)格尺寸下流固界面溫度分布Fig.8.Temperature distribution on the fluid-solid interface under different mesh sizes.

從圖7 可以看出,溫度為300 K的流體流過(guò)冷卻通道后,可以將處于900 K 熱環(huán)境的結(jié)構(gòu)外表面溫度降低到400 K 以下,充分說(shuō)明了高溫結(jié)構(gòu)中采用主動(dòng)冷卻管道的冷卻效果.

圖7 不同網(wǎng)格尺寸下管道外表面溫度分布Fig.7.Temperature distribution on the channel outer surface under different mesh sizes.

5.1.2 使用比例間距網(wǎng)格計(jì)算

為了考核FLM 在非均勻網(wǎng)格中的表現(xiàn),在上述FLM-fine1的網(wǎng)格中,允許流體域與固體域的計(jì)算點(diǎn)從各域邊界到其中線按比例變化,相鄰間距比例因子在水平與垂直方向分別為1.2 和1.5.圖9展示了流固域左端與中間部分的FLM 線系連成的局部網(wǎng)格圖.圖10 和圖11 分別為不同流速下管道外表面和流固界面的溫度沿x方向的分布曲線.為了比較,也繪出了用同等節(jié)點(diǎn)數(shù)但等間隔分布的網(wǎng)格的計(jì)算結(jié)果(FLM-uniform).

圖9 管道結(jié)構(gòu)左端與中間局部網(wǎng)格放大圖Fig.9.Enhanced meshes of the left and middle parts of the channel structure.

可以看出,用比例間距的計(jì)算結(jié)果 (FLM-ratio)與Fluent的計(jì)算結(jié)果吻合的很好,而用等間隔網(wǎng)格的計(jì)算結(jié)果(FLM-uniform)在如此少的節(jié)點(diǎn)數(shù)下的計(jì)算結(jié)果不可接受.可見(jiàn),在FLM 計(jì)算中,采用不等間距網(wǎng)格可以大量減少所需的節(jié)點(diǎn)數(shù).不過(guò),從圖10 和圖11 中的v=0.1的結(jié)果可以看出,對(duì)于流速慢的情況,即使用等間距網(wǎng)格也可以獲得滿意的計(jì)算結(jié)果.

圖10 不同流速下管道外表面溫度分布Fig.10.Temperature distribution on channel outer surface under different velocities.

圖11 不同流速下流固界面溫度分布Fig.11.Temperature distribution on fluid-solid interface under different velocities.

5.2 多根管道結(jié)構(gòu)流動(dòng)傳熱分析

本文的第2 個(gè)算例是含有3 根管道的冷卻結(jié)構(gòu),如圖12 所示.三根管道的直徑均為2 mm,其中,中間管道具有5.7°的傾斜角.整個(gè)冷卻結(jié)構(gòu)的左、右邊界除了流體的入口溫度Tin部分外其他部分均為絕熱邊界條件(q=0).三根管道中的流體的密度與比熱與算例5.1 中的一樣,均為ρ=1000 kg/m3與cp=4200 J·(kg·K)—1.計(jì)算中,整個(gè)冷卻結(jié)構(gòu)分為7 個(gè)計(jì)算域,其中第2,4,6為流體域,其他為固體域.幾何離散中,所有域沿水平方向(x方向)劃分為等間隔的101 個(gè)點(diǎn),沿垂直方向劃分為從各域壁面到其中線間距比例因子為1.2的31 個(gè)點(diǎn).圖13為FLM 分析中由所有線系連成的網(wǎng)格圖,注意到由于中間管道的傾斜使得計(jì)算域3—5的網(wǎng)格都成了非正交網(wǎng)格.

圖12 三根管道結(jié)構(gòu)尺寸與邊界條件Fig.12.Dimensions and B.C.of a three-channel structure.

圖13 三管道結(jié)構(gòu)FLM 分析線系連成的網(wǎng)格圖Fig.13.Mesh connected by line sets of all points for FLM analysis of the three channel structure.

5.2.1 流、固域采用相同條件的計(jì)算分析

首先考察整個(gè)管道結(jié)構(gòu)處于與算例5.1 所示的外部環(huán)境與邊界條件,即上下外表面所處的環(huán)境溫度為T(mén)∞=900 K,換熱系數(shù)為h=1100 W·(m2·K)—1,流體與固體的導(dǎo)熱系數(shù)λ 以及流體域的入口溫度Tin與算例5.1 中的相同.圖14為在不同流速下整個(gè)冷卻系統(tǒng)的溫度云圖,圖15為上部外表面溫度變化曲線,而圖16為斜管流體域下表面的溫度變化曲線.

從圖14 可以看出,不同流速下冷卻系統(tǒng)中間部分的溫度變化很大.當(dāng)流速很低時(shí)(如v=0.002 m/s的情況),只有流體的入口部分能得到冷卻;隨著流速的增加,冷卻效果越來(lái)越明顯;當(dāng)流速為v=0.5 m/s 時(shí),整個(gè)中間結(jié)構(gòu)都能得到冷卻.另外,從圖16 和圖14(d)可以看出,當(dāng)流速為v=0.5 m/s時(shí),斜管流體域下表面的溫度基本為流體的入口溫度,說(shuō)明在此流速下不需要中間的冷卻管道對(duì)結(jié)構(gòu)進(jìn)行冷卻.

圖14 不同流速下的溫度云圖 (a) v=0.002 m/s;(b) v=0.005 m/s;(c) v=0.01 m/s;(d) v=0.5 m/sFig.14.Contours under different velocities: (a) v=0.002 m/s;(b) v=0.005 m/s;(c) v=0.01 m/s;(d) v=0.5 m/s.

圖15 不同流速下結(jié)構(gòu)上部外表面溫度變化曲線Fig.15.Temperature variation curve on the outer surface of upper structure under different velocities.

圖16 不同流速下斜管流體域下表面溫度變化曲線Fig.16.Temperature variation curve on the lower surface of fluid domain of the oblique channel under different velocities.

5.2.2 流、固域采用不同條件的計(jì)算分析

接下來(lái)將整個(gè)結(jié)構(gòu)的上下表面對(duì)流換熱邊界條件分別設(shè)置為T(mén)∞=1000 K,h=800 W·(m2·K)—1和T∞=700 K,h=1100 W·(m2·K)—1;各計(jì)算域的導(dǎo)熱系數(shù)λ、流體域的入口溫度Tin和流速v如表1所示.圖17為計(jì)算得到的整個(gè)冷卻系統(tǒng)的溫度云圖,圖18 給出了冷卻系統(tǒng)上表面(用Upper-surface標(biāo)記)、下表面(Lower-surface 標(biāo)記)、和各區(qū)域界面上的溫度沿x方向的變化曲線,其中Interface 45 表示計(jì)算域4 與5的界面,其他類推.

表1 各計(jì)算域的材料參數(shù)與入口邊界條件Table 1.Material parameters and boundary conditions of each computational domain.

從圖17 可以看出,冷卻通道附近的溫度與冷卻流體的溫度接近,說(shuō)明冷卻流體帶走的熱量明顯.另外,從圖18 可以看出,入口(x=0)處上下外表面的溫度與環(huán)境溫度具有較大的差別,這說(shuō)明冷卻流體可以顯著地降低結(jié)構(gòu)的表面溫度,這是高溫裝備通常使用主動(dòng)冷卻熱防護(hù)的主要原因.

圖17 計(jì)算得到的管道冷卻系統(tǒng)溫度云圖Fig.17.Contours of computed temperature over the channel cooling system.

圖18 計(jì)算得到的冷卻系統(tǒng)上下表面和各區(qū)域界面的溫度變化曲線Fig.18.Variation curve of the computed temperature on the upper and lower surfaces as well as on the interfaces of the cooling system.

5.3 三維三管道結(jié)構(gòu)流固域間耦合傳熱分析

本文分析的第3 個(gè)算例是含有3 根管道的三維冷卻結(jié)構(gòu),如圖19 所示.管道的尺寸均為1.333 mm×0.667 mm.上下表面的溫度分別給定為T(mén)1=600 K 和T2=350 K.計(jì)算中,整個(gè)冷卻結(jié)構(gòu)分為6 個(gè)計(jì)算域,其中Ω4—Ω6為流體域,其他為固體域;冷卻通道內(nèi)冷卻液的入口溫度Tin和流速v見(jiàn)表2,流體的密度與比熱均為ρ=1000 kg/m3與cp=4200 J·(kg·K)—1.

圖19 含三根管道的三維冷卻結(jié)構(gòu)尺寸與邊界條件Fig.19.Dimensions and B.C.of a 3D cooling structure with three channels.

表2 三維流固結(jié)構(gòu)各計(jì)算域的材料參數(shù)與邊界條件Table 2.Material parameters and boundary conditions of each computational domain.

圖20為FLM 計(jì)算用的網(wǎng)格圖,總共有172789個(gè)配置點(diǎn),其中沿管道走向(x方向)劃分了31 個(gè)點(diǎn),靠近壁面的網(wǎng)格間距比例因子為1.2;各域沿垂直方向(z方向)均劃分為11 個(gè)點(diǎn),壁面比例因子為1.5;沿橫向(y方向),除上下蓋板的外伸部分為25 個(gè)點(diǎn)外,其他部分均為21 個(gè)點(diǎn),壁面比例因子均為1.5.圖21為橫截面網(wǎng)格與局部放大圖,可以看出下蓋板的外伸部分是曲面,而且靠近右下角的網(wǎng)格極度非正交.

圖20 所有配置點(diǎn)的線系連成的網(wǎng)格Fig.20.Mesh connected by line sets of all collocation points.

圖21 橫截面上的配置點(diǎn)線系連成的網(wǎng)格與局部放大圖Fig.21.Mesh connected by line sets on the transverse section and a locally refined part around a corner.

圖22 是計(jì)算得到的總體結(jié)構(gòu)上的溫度云圖,而圖23 是通過(guò)三根通道中線的垂直面上的溫度云圖.圖24 和圖25 分別為沿圖22 所示的線段L1與L2上的溫度分布.為了比較,也用商業(yè)有限容積法軟件FLUENT 和有限單元法軟件COMSOL 對(duì)該問(wèn)題進(jìn)行了計(jì)算,使用的網(wǎng)格是達(dá)到網(wǎng)格收斂下的網(wǎng)格.對(duì)比可以看出,FLM 計(jì)算結(jié)果與其他兩種方法的結(jié)果吻合的很好,證明了本文方法在使用非規(guī)則網(wǎng)格方面的正確性,這是傳統(tǒng)的有限差分法難以實(shí)現(xiàn)的.另外,表3 列出了采用3 種方法計(jì)算該問(wèn)題時(shí)所用的節(jié)點(diǎn)數(shù)與時(shí)間的比較.其中,COMSOL只列了一種網(wǎng)格的結(jié)果,因?yàn)榫W(wǎng)格稍微密一點(diǎn)時(shí)計(jì)算不收斂,試算發(fā)現(xiàn)這主要是因?yàn)榱黧w與固體的導(dǎo)熱系數(shù)相差較大所致.從表3 可以看出,在同等精度下FLM 需用的網(wǎng)格數(shù)最少,消耗時(shí)間也較少;COMSOL 用的時(shí)間最多,而且在流體入口處附近與其他兩種方法的結(jié)果有較大的差別,穩(wěn)定性也差,這也許是有限單元法在處理流體相關(guān)問(wèn)題中存在的固有弱點(diǎn)[20].此外,FLM 與FLUENT 用粗、細(xì)兩種網(wǎng)格算的結(jié)果基本一樣,說(shuō)明二者有很好的穩(wěn)定性.

圖22 總體結(jié)構(gòu)上的溫度云圖Fig.22.Contour plot of computed temperature over the entire structure.

圖23 管道走向垂直中面上的溫度云圖Fig.23.Contour plot of computed temperature over the vertical middle plan.

圖24 沿圖22 中所示的線段 L1 上的溫度分布Fig.24.Computed temperature along line L1 marked in Fig.22.

圖25 沿圖22 中所示的線段 L2 上的溫度分布Fig.25.Computed temperature along line L2 marked in Fig.22.

表3 三種計(jì)算方法的節(jié)點(diǎn)數(shù)與計(jì)算時(shí)間比較Table 3.Comparison of total number of nodes and computational time for three methods.

6 討論與結(jié)論

論文提出了一種全新的數(shù)值計(jì)算方法—有限線法,并將其應(yīng)用于求解流體域與固體域之間的耦合傳熱問(wèn)題,根據(jù)算法的理論與在冷卻通道結(jié)構(gòu)中的算例分析,可以得出如下的結(jié)論:

1)本文通過(guò)線段方向?qū)?shù)建立的求導(dǎo)公式可以用統(tǒng)一的格式計(jì)算多維問(wèn)題的任意高階的空間偏導(dǎo)數(shù);

2)創(chuàng)建的有限線法是一種配點(diǎn)型計(jì)算方法,其用法與傳統(tǒng)的有限差分法一樣方便、靈活,而且能模擬曲的和不規(guī)則的幾何形狀;

3)有限線法可以直接由問(wèn)題的控制微分方程和邊界條件建立系統(tǒng)方程組,不需要變分原理、能量原理及其他需要積分的數(shù)學(xué)操作;

4)本文雖然以流固耦合傳熱問(wèn)題為研究背景,但導(dǎo)出的公式也適用于求解其他工程問(wèn)題的偏微分方程邊值問(wèn)題,是一種通用的數(shù)值計(jì)算方法;

5)因?yàn)橛邢蘧€法是一種強(qiáng)形式的計(jì)算方法,通常需要采用高階線系計(jì)算,經(jīng)驗(yàn)表明,采用四階精度的中心線系和二階精度的迎風(fēng)線系可以獲得高精度和高效率的計(jì)算結(jié)果;

本文所述方法是關(guān)于空間偏導(dǎo)數(shù)的離散的方法,考慮的控制方程是不含時(shí)間項(xiàng)的定常問(wèn)題,但所導(dǎo)出的公式也同樣適用于含有時(shí)間項(xiàng)的非定常問(wèn)題,此時(shí),時(shí)間項(xiàng)可采用成熟的方法進(jìn)行離散[24].